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Supercritical  Combustion  of  Liquid  Oxygen  and  Hydrocarbon  for 
Staged-Combustion  Cycle  Engine  Technology  Development 

SUMMARY 

An  integrated  modeling  and  numerical  program  has  been  conducted  to 
substantially  improve  the  fundamental  knowledge  of  supercritical  fluid  dynamics  and 
combustion  of  oxygen  and  hydrocarbon  fuels  under  conditions  representative  of 
contemporary  liquid-propellant  rocket  engines.  Emphasis  was  placed  on  the  swirler 
injector  flow  dynamics  and  flame  stabilization  of  LOX  and  methane.  The  fundamental 
characteristics  of  counter-flow  diffusion  flames  of  oxygen  and  hydrogen  (or  methane) 
were  also  explored.  The  analysis  was  based  on  the  complete  conservation  equations  for  a 
multi-component  chemically  reacting  mixture,  and  accommodated  general-fluid 
thermodynamics  and  transport  phenomena  valid  for  the  entire  range  of  fluid  states  of 
concern.  Turbulence  closure  was  treated  using  the  large-eddy-simulation  (LES) 
technique.  Start-of-the-art  closure  schemes  for  subgrid-scale  dynamics  and 
turbulence/chemistry  interactions  were  implemented. 

The  work  addressed  the  following  issues: 

1.  cryogenic  fluid  dynamics  of  pressure  swirl  injectors  at  supercritical 
conditions, 

2.  counterflow  diffusion  flames  of  general  fluids,  including  oxygen/hydrogen 
and  oxygen/methane  mixtures,  and 

3.  flame  stabilization  and  spreading  of  liquid  oxygen  (LOX)  and  methane  at 
supercritical  conditions. 

Various  underlying  physiochemical  mechanisms  associated  with  swirl  co-axial 
injectors  have  been  studied  in  detail.  These  include  flow  evolution,  flame  stabilization 
and  spreading,  heat  transfer,  and  acoustic  response.  The  effects  of  design  attributes  (e.g., 
injection  port  size  and  location,  center  post  recess  distance,  etc.)  and  operating  conditions 
(e.g.,  chamber  pressure,  velocity,  and  temperature,  swirl  strength,  etc)  on  injector 
characteristics  were  assessed.  Results  not  only  enhance  the  basic  understanding  of  the 
subject  problems,  but  also  provide  a  quantitative  basis  to  identify  and  prioritize  the  key 
design  parameters  and  flow  variables  that  exert  dominant  influence  on  the  injector 
behavior  in  different  environments. 

Results  from  this  research  project  have  led  to  the  following  publications. 

1.  “Cryogenic  Fluid  Dynamics  of  Pressure  Swirl  Injectors  at  Supercritical 
.  Conditions,”  by  N.  Zong  and  V.  Yang,  Physics  of  Fluids ,  Vol.  20,  2008, 
056103  (1-14). 
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2.  “Counterflow  Diffusion  Flames  of  General  Fluids:  Oxygen/Hydrogen 
Mixtures,”  by  G.  Ribert,  N.  Zong,  V.  Yang,  L.  Pons,  N.  Darabiha,  and  S. 
Candel,  Combustion  and  Flame,  Vol.  154,  2008,  pp.  319-330. 

3.  “Mass  Transfer  and  Combustion  in  Transcritical  Non-Premixed 
Countcrflows,”  by  L.  Pons,  N.  Darabiha,  S.  Candel,  G.  Ribert,  and  V.  Yang, 
Combustion  Theory  and  Modeling,  Vol.  13,  No.  1, 2009,  pp.  57-81. 

4.  “Supercritical  Combustion  of  LOX/Methane  Stabilized  by  a  Splitter  Plate,”  by 
N.  Zong,  G.  Ribert,  and  V.  Yang,  AIAA  Paper  2007-0575,  presented  at  the 
45th  AIAA  Aerospace  Sciences  Meeting  and  Exhibit,  January  2007. 

5.  “Large  Eddy  Simulation  of  Combustion  of  Liquid  Oxygen  and  Methane  in  a 
Supercritical  Environment,”  by  G.  Ribert,  N.  Zong,  and  V.  Yang,  presented  at 
European  Community  Conference  on  Large  Eddy  Simulation  for  Advanced 
Design  of  Combustion  Systems,  Rouen,  France,  May  2007. 

6.  “A  Flamelet  Approach  for  Modeling  of  Liquid  Oxygen  (LOX)/Methane 
Flames  at  Supercritical  Pressures,”  by  N.  Zong,  G.  Ribert,  and  V.  Yang, 
AIAA  Paper  2008-946,  presented  at  46th  AIAA  Aerospace  Sciences  Meeting 
and  Exhibit,  January  2008. 

7.  “Cryogenic  Fluid  Dynamic  Response  of  Swirl  Injector  to  External  Forcing  at 
Supercritical  Conditions,”  by  H.  Huo,  N.  Zong,  and  V.  Yang,  AIAA  Paper 
2009-0233,  presented  at  47th  Aerospace  Sciences  Meeting,  January  2009. 
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Task  1 


Cryogenic  Fluid  Dynamics  of  Pressure  Swirl  Injectors  at  Supercritical  Conditions 


Summary 


A  comprehensive  numerical  analysis  has  been  conducted  to  explore  the 
development  of  liquid  oxygen  flow  in  pressure  swirl  injectors  operating  at  supercritical 
pressures.  The  model  is  based  on  full-conservation  laws  and  accommodates  real-fluid 
thermodynamics  and  transport  phenomena  over  the  entire  range  of  fluid  states  of  concern. 
Three  different  flow  regimes  with  distinct  characteristics,  the  developing,  stationary,  and 
accelerating  regimes,  are  identified  within  the  injector.  Results  are  compared  with 
predictions  from  classical  hydrodynamic  theories  to  acquire  direct  insight  into  the  flow 
physics  involved.  In  addition,  various  flow  dynamics  arc  investigated  by  means  of  the 
spectral  and  proper-orthogonal-decomposition  techniques.  The  interactions  between  the 
hydrodynamic  instabilities  in  the  LOX  film  and  acoustic  oscillations  in  the  gaseous  core 
are  clearly  observed  and  studied.  The  influences  of  flow  conditions  (mass  flowrate,  swirl 
strength  of  the  injected  fluid,  and  ambient  pressure)  and  injector  geometry  (injector 
length  and  tangential  entry  location)  on  the  injector  flow  behavior  are  characterized 
systematically  in  terms  of  the  LOX  film  thickness  and  spreading  angle.  The  axial  and 
azimuthal  momentum  exchange  and  loss  mechanisms  are  also  examined. 

1.1  Introduction 

This  paper  deals  with  the  cryogenic  fluid  dynamics  of  a  pressure  swirl  injector,  as  shown 
schematically  by  the  center  element  of  a  co-axial  injector  in  Fig.  1.  The  configuration  is 
representative  of  contemporary  liquid-propellant  rocket  injectors  for  liquid/liquid  and 
gas/liquid  mixtures.1  Liquid  propellant  is  introduced  tangentially  into  the  center  post,  and 
then  forms  a  swirling  film  attached  to  the  wall  due  to  centrifugal  force.  A  hollow  gas 
core  exists  in  the  center  region  in  accordance  with  the  conservation  of  angular 
momentum.  The  film  exits  the  injector  as  a  thin  sheet  and  impinges  onto  the  surrounding 


breakup  length 


fuel 


liquid  oxygen 


Fig.  1  Schematic  of  swirl  co-axial  injector. 
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fuel  stream.  The  swirl  injection/atomization  process  basically  involves  two  mechanisms: 
disintegration  of  the  liquid  sheet  as  it  swirls  and  stretches;  and  sheet  break-up  due  to  the 
interaction  with  the  surrounding  co-axial  flow.  In  spite  of  the  broad  applications  of  swirl 
co-axial  injectors,  fundamental  studies  on  the  mixing  and  combustion  processes  appear  to 
be  limited,  especially  at  scales  sufficient  to  identify  the  underlying  mechanisms  dictating 
the  flow  evolution  and  flame  dynamics,  and  at  high-pressures. 

Compared  with  jet  injectors,  swirl  injectors  in  liquid-propellant  rocket  engines 
(LREs)  distinguish  themselves  in  several  aspects.  First,  the  non-uniform  mixing  of 
propellants  in  the  jet  core  region  is  avoided  and  the  intra-clemcnt  mixing  is  significantly 
improved  because  of  the  outward  spreading  of  the  liquid  spray.  High  mixing  efficiency 
is  thus  possible  even  for  a  large  injector  flowrate.  Second,  the  large  flow  passage  in  a 
swirl  injector  renders  the  atomization  characteristics  less  sensitive  to  manufacturing 
errors.  The  injector  is  also  less  susceptible  to  choking  and  cavitation.  Third,  the  injected 
fluid  is  discharged  into  the  chamber  as  a  hollow  spray  cone.  The  thickness  of  the  liquid 
film  becomes  thinner  as  it  swirls  and  spreads  outward.  The  resultant  mean  diameter  of 
droplets  is  2.2  to  2.5-times  smaller  than  that  produced  by  a  jet  injector  with  the  same 
pressure  drop  and  mass  flowrate.  The  droplet  size  distribution  is  also  more  uniform. 
Finally,  rocket  swirl  injectors  feature  large  aspect  ratios  of  up  to  20,  due  to  the 
manifolding  considerations  of  propellant  supply.  The  viscous  loss  along  the  injector  wall 
exerts  significant  influence  on  the  flow  evolution  and  subsequently  alters  liquid 
atomization  characteristics. 

Hulka  et  al.  4  conducted  a  series  of  cold-flow  studies  to  optimize  the  design  of 
swirl  co-axial  injectors  for  liquid  rockets.  Three  different  types  of  experiments  were 
performed  in  a  quiescent,  atmospheric  environment.  Still  photographs  were  first  taken  to 
measure  the  spreading  angle  of  a  water  spray  without  a  surrounding  co-flow.  The 
circumferential  spray  uniformity  and  droplet  size  distribution  were  then  measured  under 
conditions  with  and  without  co-flow  nitrogen  gas  by  means  of  a  grid  pattemator  and  a 
Malvern  droplet  size  analyzer,  respectively.  The  co-flow  reduces  the  mean  droplet  size, 
but  widens  the  droplet  size  distribution.  Measurements  were  also  made  of  the  Rupc 
mixing  efficiency  using  water  and  a  sucrose  solution  to  simulate  liquid  oxygen  and 
gaseous  hydrogen,  respectively.  A  broad  range  of  fuel-to-oxidizer  mixture  (0.94-17.8) 
and  velocity  (1.15-4.28)  ratios  at  one  atmosphere  was  considered.  The  mixing  efficiency 
could  be  greatly  enhanced  by  increasing  the  initial  swirl  strength  of  the  oxidizer,  which  is 
a  function  of  the  injector  geometric  characteristic  parameter.  The  oxidizer  mass  flowTatc 
has  little  influence  on  the  mixing  efficiency.  Thus,  an  oxidizer  swirl  element  with  a  mass 
flowrate  greater  than  that  of  a  shear  co-axial  injector  can  still  achieve  better  mixing 
efficiency. 

A  similar  cold-flow  experiment  with  water  and  nitrogen  as  injectants  was 
performed  by  Sasaki  et  al.5  at  room  conditions.  Special  attention  was  given  to  the  effect 
of  the  center  post  recess,  which  tends  to  narrow  the  spreading  angle  and  cause  a 
deformation  of  the  spray  cone.  The  influence  of  the  annular  co-flow  on  spray 
characteristics  was  also  investigated.  The  swirling  liquid  sheet  that  blocked  the  outer 
gas-flow  passage  in  a  recessed  injector  was  blown  off  in  a  mushroom  shape  accompanied 
with  a  high-intensity  screaming  sound.  This  phenomenon,  known  as  self-pulsation,  may 
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result  in  reduced  combustion  efficiency  and  high-amplitude  pressure  oscillations  in  the 
chamber. 1  0 

The  recess  effects  of  the  center  element  on  the  mixing  characteristics  of  a  swirl 
co-axial  injector  were  also  examined  by  Han  et  al.  and  Kim  et  al.  using  kerosene  as  the 
fuel  and  water  as  the  oxidizer  simulant.  Backlit  stroboscopic  photographs  were  taken  to 
determine  the  liquid  spreading  angle  and  breakup  length.  The  median  droplet  size, 
propellant  mass  distribution,  and  mixing  efficiency  were  measured  using  a  phase  Doppler 
particle  analyzer  (PDPA)  and  a  mechanical  pattemator.  The  influence  of  the  inner- 
stream  spreading  angle  was  taken  into  account  by  considering  four  different  recess 
numbers  in  the  range  of  0.71-1.37,  defined  as  the  ratio  of  the  recess  length  to  the  distance 
from  the  center  post  to  the  position  where  the  swirling  liquid  sheet  impinges  onto  the 
outer  wall  of  the  annular  passage.  The  mixing  efficiency  and  propellant  mass  distribution 
were  found  to  be  very  sensitive  to  the  recess  length.  The  Sauter  mean  diameter  of 
droplets  decreased  slightly  with  increasing  recess  length,  and  could  be  correlated  with  the 
empirical  equations  suggested  by  Lefebvre9  and  Jones.10 

Marchione  et  al."  recently  used  a  PDPA  system  to  measure  the  mean  droplet  size 
and  velocity  profiles  of  a  kerosene  spray  produced  by  a  swirl  injector.  The  experiment 
was  conducted  in  a  quiescent  environment  at  one  atmosphere.  Instantaneous  images 
were  taken  with  a  high-speed  CCD  camera  to  analyze  both  the  static  and  oscillatory 
behaviors  of  the  spray.  The  spray  angle  estimated  based  on  the  measured  mean  droplet 
size  and  velocity  agreed  well  with  the  results  derived  from  the  images.  Two  oscillation 
modes  of  the  spray  at  a  low  frequency  of  100  Hz  and  a  high  frequency  of  1800  Hz  were 
identified  by  analyzing  the  Mie  scattering  intensity  in  different  regions.  Those  modes  of 
oscillation  could  induce  a  variation  of  the  local  air/fuel  ratio. 

All  the  aforementioned  studies2  5’  ,8,11  were  conducted  at  one  atmosphere  without 
considering  the  effects  of  the  elevated  pressures  typically  encountered  in  operational 
liquid-propellant  rocket  engines.  Strakey  et  al.12  investigated  the  spray  characteristics  of 
swirl  co-axial  injectors  over  a  broad  range  of  liquid-to-gas  momentum  ratios  (0.1-100)  at 
a  chamber  pressure  of  2.97  MPa  using  water  and  helium/nitrogen  as  propellant  simulants. 
The  spreading  angle  of  the  liquid  spray  decreases  with  the  increase  of  the  co-axial  gas 
momentum,  and  appears  almost  identical  to  that  produced  by  a  shear  co-axial  injector  at 
the  lowest  liquid-to-gas  momentum  ratio  (0.1).  The  mean  droplet  size,  however, 
becomes  smaller  compared  with  that  of  a  shear  co-axial  injector  at  the  same  flowrate. 
Recently,  Yoon  and  coworkers  studied  the  effects  of  ambient  gas  density  on  the  liquid- 
sheet  spreading  angle  and  breakup  length  of  single  swirl"  and  swirl  co-axial1 4  injectors. 
Water  was  employed  as  the  oxidizer  simulant.  The  measured  spray  angle  prior  to  the 
sheet  breakup  remains  almost  fixed  in  the  pressure  range  of  1  -40  atm.  The  sheet  breakup 
length,  however,  decreases  with  increasing  pressure  due  to  the  enhanced  aerodynamic 
force  as  the  ambient  gas  density  increases.  Hla 15  conducted  hot-fire  experiments  for 
multi-element  swirl  injectors  with  liquid  oxygen  and  gaseous  hydrogen  at  mixture  ratios 
of  5. 2-6. 9  and  a  chamber  pressure  of  10.3  MPa.  Results  indicated  that  the  C*  efficiency 
increases  with  increasing  mixture  ratio.  Hot-fire  studies  of  liquid  oxygen  and  gaseous 
hydrogen  were  also  performed  by  Tamura  and  colleagues10,1  for  both  single-  and  multi¬ 
element  swirl  injectors  at  moderate  chamber  pressures  (2. 6-3. 5  MPa)  and  mixture  ratios 
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of  4.0-8. 0.  It  was  reported  that  an  enhanced  C*  efficiency  could  be  obtained  at  a  greater 
mixture  ratio. 

In  addition  to  experimental  studies,  numerical  simulations  were  recently 
performed  to  explore  the  flow  evolution  of  swirl  co-axial  injectors.  Kim  and  Heistcr1" 
investigated  hydrodynamic  instabilities  of  the  swirling  liquid  jet  within  the  recessed 
region  of  a  swirl  co-axial  injector.  A  locally  homogeneous  flow  model,  previously 
developed  for  treating  flow  instabilities  in  shear  co-axial  injectors,19  was  adopted,  by 
assuming  the  thermodynamic  and  momentum  equilibria  between  the  gas  and  liquid 
phases  throughout  the  entire  flowfield.  The  approach  also  employed  an  incompressible- 
flow  assumption  with  pre-specified  gas  and  liquid  densities.  Real-fluid  thermodynamics 
and  property  variations  inherent  in  operational  injector  flows  were  not  taken  into  account, 
to  simplify  the  analysis.  Initial  results  indicated  that  the  swirling  motion  and  fluctuation 
of  the  injected  liquid  became  more  pronounced  with  an  increase  in  the  gas-to-liquid 
density  ratio.  Park  and  Heister2"  simulated  the  free  surface  and  spray  shape  of  injected 
liquid  in  a  swirl  injector.  The  analysis  was  based  on  a  boundary  element  method  along 
with  the  assumption  of  an  aixsymmetrical,  inviscid,  and  incompressible  flow.  Calculated 
film  thickness  and  spray  cone  angle  agreed  well  with  those  predicted  by  the  classical 
theory.  Canino  et  al.  1  studied  the  flow  development  behind  the  center  post  in  a  swirl  co¬ 
axial  injector.  Emphasis  was  placed  on  the  vortex-shedding  characteristics  due  to  the 
bluff-body  dynamics  within  the  injector  passage. 

Although  much  useful  information  has  been  obtained  about  the  liquid  atomization 
and  mixing  characteristics  of  swirl  co-axial  injectors,  most  existing  studies  were 
conducted  either  at  room  temperatures  using  water  as  the  oxidizer  simultant,2  5,7,8,11  14  or 
at  hot- fire  conditions1 5-17  where  the  available  diagnostics  are  limited  by  the  difficulties  in 
probing  the  flowfield  in  such  a  high-pressure,  high-temperature  environment.  Very 
limited  effort  has  been  applied  to  investigate  detailed  injector  flow  dynamics  under 
conditions  representative  of  operational  rocket  engines,  where  the  chamber  pressure  often 
exceeds  the  thermodynamic  critical  pressure  of  the  injected  fluid.  Furthermore,  the  flow 
evolution  within  the  injector  has  not  been  carefully  explored,  except  for  the  analytic 
studies  based  on  classical  hydrodynamics  theories  discussed  in  Refs.  1  and  and 

1 X 

numerical  studies  based  on  an  incompressible-flow  assumption. 

The  purpose  of  the  present  work  is  to  remedy  these  deficiencies  by  conducting 
high-fidelity  simulations  of  cryogenic  fluid  dynamics  of  swirl  injectors  at  supercritical 
pressures,  mimicking  the  configurations  and  flow  conditions  of  contemporary  rocket 
engines  using  liquid  oxygen  (LOX)  as  the  oxidizer.'4  The  formulation  accommodates 
full  conservation  laws,  and  takes  into  account  real-fluid  thermodynamics  and  transport 
phenomena.  The  specific  objectives  of  the  study  are:  1)  to  characterize  the  detailed  flow 
physics  in  the  injector  flow  path;  2)  to  explore  various  underlying  mechanisms  dictating 
the  fluid  atomization  and  energy-transfer  behaviors;  and  3)  to  identify  and  prioritize  key 
injector  design  attributes  and  operating  conditions  critical  to  the  injector  performance. 

1.2  Theoretical  Formulation 

The  basis  of  the  present  study  is  the  general  theoretical  framework  for  treating 
supercritical  fluid  injection  and  mixing,  as  detailed  by  Meng  et  al.24  and  Zong  et  al. 
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The  approach  has  been  used  to  investigate  the  transport  and  dynamics  of  droplets,”4  jet 
mixing,25  and  chemically  reacting  shear  flows2'1  over  a  broad  range  of  thermodynamic 
states  and  flow  conditions.  In  brief,  the  formulation  accommodates  complete 
conservation  equations  of  mass,  momentum,  energy,  and  species  transport.  Full  account 
is  taken  of  general-fluid  thermodynamics  and  transport  phenomena  over  the  entire 
temperature  and  pressure  regimes  of  concern.  Turbulence  closure  is  achieved  by  means 
of  a  large-eddy-simulation  (LES)  technique,  in  which  large-scale  motions  are  calculated 
explicitly,  and  the  effects  of  unresolved  small-scale  turbulence  are  modeled  either 
analytically  or  empirically.  The  governing  Favre-filtered  conservation  equations  are 
derived  by  filtering  small-scale  dynamics  from  resolved  scales  over  a  well-defined  set  of 
spatial  and  temporal  intervals.2  ’  The  effects  of  subgrid-scale  (sgs)  motions  are  treated 
using  the  approach  proposed  by  Erlebacher  ct  al,29  which  employs  a  Favre-averaged 
generalization  of  the  Smagorinsky  eddy  viscosity  model.  The  Smagorinsky  coefficient 
CR  («0.01)  and  C,  (~  0.007)  are  determined  empirically.”  It  should  be  noted  that 
accurate  modeling  of  sgs  dynamics  under  supercritical  conditions  remains  a  challenging 
task,  due  to  complications  arising  from  rapid  property  variations  and  real-fluid 
thermodynamics  and  transport.  Very  limited  effort  has  been  applied  so  far  to  quantify  the 
effects  of  real-fluid  thermodynamics  on  small-scale  turbulent  structures  and  a  well- 
calibrated  sgs  model  for  supercritical  fluid  flows  is  currently  not  available.  This  issue 
will  be  addressed  in  our  subsequent  work. 

The  fluid  volumetric  behavior  is  evaluated  using  a  modified  Soave-Redlich- 
Kwong  equation  of  state11  because  of  its  validity  over  a  broad  range  of  fluid  states  and 
simplicity  of  numerical  implementation.  Thermodynamic  properties,  such  as  enthalpy, 
Gibbs  energy,  and  constant-pressure  specific  heat,  are  derived  directly  from  fundamental 
thermodynamic  theories.  They  are  expressed  as  the  sum  of  an  ideal-gas  property  at  the 
same  temperature  and  a  thermodynamic  departure  function  accounting  for  dense-fluid 
correction. 3”  Transport  properties,  such  as  viscosity  and  thermal  conductivity,  are 
estimated  using  an  extended  corresponding-state  theory 33  14  along  with  a  32-term 
Benedict-Webb-Robin  equation  of  state.  Mass  diffusivity  is  obtained  by  the  Takahashi 
method  calibrated  for  high  pressure  conditions. 3<  The  implementation  and  validation  of 
the  property  evaluation  schemes  are  detailed  in  Refs.  24  and  32. 

1.3  Numerical  Framework 

The  theoretical  formulation  outlined  above  is  numerically  solved  by  means  of  a 
preconditioning  scheme  incorporating  a  unified  treatment  of  general  fluid 
thermodynamics.3'1  All  the  numerical  properties,  including  the  preconditioning  matrix, 
Jacobian  matrices,  and  eigenvalues,  are  derived  directly  from  fundamental 
thermodynamics  theories,  rendering  a  self-consistent  and  robust  algorithm.  The 
numerical  formulation  can  accommodate  any  equation  of  state,  and  is  valid  for  fluid 
flows  at  all  speeds  and  at  all  fluid  thermodynamic  states.  Further  efficiency  is  achieved 
by  employing  temperature  instead  of  enthalpy  as  the  primary  dependent  variable  in  the 
preconditioned  energy  equation. 3  This  procedure  eliminates  laborious  iterations  in 
solving  the  equation  of  state  for  temperature,  and  consequently  facilitates  the  load 
balance  among  computational  blocks  in  a  distributed  computing  environment.  The 
resultant  scheme  is  highly  efficient  and  suitable  for  parallel  processing. 
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The  numerical  framework  employs  a  density-based,  finite-volume  methodology 
along  with  a  dual-time-step  integration  technique. 38  Temporal  discretization  is  obtained 
using  a  second-order  backward  differencing  scheme,  and  the  inner-loop  pseudo-time  term 
is  integrated  with  a  four-step  Runge-Kutta  scheme.  Spatial  discretization  is  achieved 
with  a  fourth-order,  central-difference  scheme  in  generalized  coordinates/0  A  nine-point 
stencil  is  employed  to  evaluate  the  convective  flux  in  each  spatial  direction  to  improve 
the  spectral  resolution  of  small-scale  turbulence  structures.  Fourth-order  scalar 
dissipation  with  the  coefficient  e4  =0.001  is  applied  to  ensure  numerical  stability  with 
minimum  contamination  of  the  solution. 

The  overall  accuracy  of  the  present  scheme  within  the  context  of  LES  was 
carefully  assessed  based  on  the  decay  of  the  kinetic  energy  of  isotropic  turbulence. 
Calculations  were  performed  for  an  isotropic  turbulent  flow  in  a  cubic  box  of  a  non- 
dimensional  width  of  27i.  The  experimental  results  of  Comet-Bellot  and  Corrsin  (CBC)40 
were  selected  as  the  benchmark  with  an  initial  Taylor  Reynolds  number  of  80  on  a 
32  x  32  x  32  grid.  The  predicted  temporal  evolution  of  turbulent  kinetic  energy  agrees 
well  with  experimental  data.  The  effects  of  numerical  and  sgs  dissipation  on  the 
evolution  of  turbulent  kinetic  energy  were  further  assessed  by  either  turning  those  terms 
on  or  off  or  reducing  the  corresponding  coefficients  by  one  half.  Results  indicated  that 
the  dissipation  associated  with  the  sgs  model  overshadows  that  of  the  numerical  scheme, 
and  the  numerical  dissipation  only  serves  to  maintain  numerical  stability. 

The  present  flowfield  involves  exceedingly  large  property  gradients  between  the 
injected  LOX  and  ambient  gaseous  oxygen.  The  density  ratio  may  reach  a  level  of  10  for 
a  chamber  pressure  100  atm.  A  second-order  scalar  dissipation  with  a  total-vanation- 
diminishing  switch  developed  by  Swanson  and  Turkel41  was  thus  applied  to  suppress 
numerical  oscillations  in  regions  with  steep  property  variations.  Care  was  also  exercised 
to  ensure  sufficient  grid  resolution  for  properly  capturing  the  hydrodynamic  instability 
waves  inherent  in  the  LOX  film. 

Finally,  a  multiblock  domain  decomposition  technique,  along  with  static  load 
balance,  is  employed  to  facilitate  the  implementation  of  parallel  computation  with 
message  passing  interfaces  (MPI)  at  the  domain  boundaries.  The  parallelization 
methodology  is  robust  and  the  speedup  is  almost  linear. 

1.4  Injector  Configuration  and  Boundary  Conditions 

The  dynamics  of  liquid  oxygen  (LOX)  injected  through  a  pressure  swirl  injector 
into  a  free  volume  preconditioned  with  gaseous  oxygen  at  a  supercritical  pressure  is 
investigated.  Figure  2  shows  the  physical  model  considered  herein.  It  consists  of  three 
major  parts:  tangential  inlets,  a  vortex  chamber,  and  a  discharge  nozzle.1  6  The  baseline 
geometry  and  mass  flowrate  are  summarized  in  Table  1,  where  R$,  Rn,  and  Rp  denote 

the  radii  of  the  vortex  chamber,  discharge  nozzle,  and  tangential  inlet,  respectively,  and 
L  the  injector  length.  The  parameters  are  chosen  to  be  identical  to  those  of  the  center 
element  of  the  swirl  co-axial  injector  in  the  RD-01 10  liquid  rocket  engine,  the  third-stage 
engine  for  the  Soyuz  space  launch  vehicle.42  The  geometric  characteristic  constant,  K , 
which  determines  the  initial  swirl  strength  of  the  injected  fluid,1  is  defined  as 
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K  s  AnRin  I  AinRn ,  with  An  being  the  discharge  nozzle  area,  Ain  the  total  area  of  the 
tangential  inlets,  and  Rin  the  radial  location  of  the  tangential  entry. 


Fig.  2  Schematic  of  pressure  swirl  injector. 
Table  1  Baseline  geometry  and  mass  flowrate 


Rs  (mm) 

Rn  (mm) 

Rp  (mm)  L  (mm) 

K 

m  (kg/'s) 

2.5 

2.5 

0.85  25 

3.2 

0.15 

The  computational  domain  includes  the  interior  of  the  swirl  injector  and  a 
downstream  region  measuring  8 Rs  and  40/?t  in  the  radial  and  axial  directions, 
respectively.  Because  of  the  enormous  computational  effort  required  for  calculating  the 
flow  evolution  in  the  entire  regime,  only  a  cylindrical  sector  with  periodic  boundary 
conditions  specified  in  the  azimuthal  direction  is  treated.  Such  an  axisymmetric 
simulation  introduces  several  limitations:  1)  the  tangential  inlets,  which  consist  of  six 
small  circular  ports  on  the  injector  wall,  are  simplified  to  a  1  mm  wide  slit  on  the  radial 
boundary  of  the  injector;  2)  the  flow  variations  in  the  azimuthal  direction  are  neglected; 
3)  the  vortex  stretching  mechanism,  which  is  responsible  for  the  energy  transfer  from 
large-  to  small-scale  structures  in  the  flowfield,  is  not  considered.  In  spite  of  those 
limitations,  previous  studies"*’’41  have  indicated  that  the  present  numerical  analysis  is  able 
to  capture  many  unique  mechanisms  dictating  supercritical  fluid  injection  and  mixing 
dynamics,  including  thermodynamic  non-idealities,  density  stratification,  interfacial 
instability,  and  baroclinic  vorticity  production.  In  addition,  various  important  unsteady 
flow  features,  such  as  the  interactions  of  hydrodynamic  instabilities  and  acoustic 
oscillations  within  the  swirl  injector,  are  well  explored  in  the  present  work.  The  results 
provide  a  quantitative  basis  for  identifying  the  key  processes  and  injector  parameters  that 
determine  the  liquid  sheet  breakup  and  atomization  characteristics  under  supercritical 
pressures. 

At  the  inlet,  the  bulk  radial  and  azimuthal  velocities  are  selected  to  match  the 
mass  flowrate,  m ,  and  the  swirl  strength  of  the  injected  fluid,  which  is  estimated  based 
on  the  injector  geometric  constant,  K.  The  fluid  temperature  is  fixed  and  pressure  is 
obtained  through  a  one-dimensional  approximation  to  the  radial  momentum  equation 
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(i.e.,  dp/dr  =  -p(dur/dt)-ur(dur/dr)  +  nl/r  ).  Turbulence  is  provided  by 
superimposing  broadband  noise  onto  the  mean  velocity  profiles.  The  disturbances  are 
produced  by  a  Gaussian  random-number  generator  with  an  intensity  of  8%  of  the  mean 
quantities,  sufficient  to  trigger  the  instabilities  inherent  in  the  flowfield.  A  white  noise 
with  a  higher  turbulence  intensity  (10%)  is  also  considered  for  the  baseline  case  (Case  1 
in  Table  2).  No  obvious  changes  in  the  injector  flow  dynamics  are  observed  in  terms  of 
the  spectral  characteristics  of  flow-oscillations  at  different  locations.  At  the  downstream 
boundary,  extrapolation  of  primitive  variables  from  the  interior  may  cause  undesired 
reflection  of  numerical  waves  propagating  into  the  computation  domain.  Thus,  the 
characteristic  boundary  conditions  proposed  by  Poinsot  and  Lele44  for  ideal  gases  and 
extended  to  real-fluid  flows17  are  incorporated  into  the  present  preconditioning  scheme, 
and  a  time-invariant  back  pressure  is  specified.  At  the  radial  boundary,  the  temperature 
and  velocity  components  are  extrapolated  from  the  interior  with  a  fixed  ambient  pressure. 
Because  the  computational  domain  is  sufficiently  large,  the  truncation  effects  of  the 
radial  and  downstream  boundaries  on  the  injector  near-field  flow  dynamics  arc 
negligible.  Finally,  the  non-slip  adiabatic  conditions  are  enforced  along  the  solid  wall. 

1.5  Results  and  Discussion 

The  liquid  film  thickness,  h,  and  spreading  angle,  2a,  at  the  injector  exit  arc 
often  employed  to  characterize  the  performance  of  a  swirl  injector.'  The  former  dictates 
the  size  of  the  fluid  parcel  after  the  film  breaks  up,  and  the  latter  affects  the  intra-element 
mixing  efficiency.  The  spreading  angle  is  defined  as  twice  the  apex  angle  of  the 
asymptotic  cone  to  the  hyperboloid  of  revolution  corresponding  to  the  profile  of  the 
spray,  and  can  be  calculated  by  the  ratio  of  the  circumferential  to  the  axial  velocities  at 
the  injector  exit,  i.e.,  a  =  tan  ' (ue0 / ua) The  key  variables  influencing  the  injector  flow 
dynamics  include  the  geometric  constant,  K ,  injector  length,  L ,  injector  exit  diameter, 
Dn ,  and  thermophysical  properties  of  the  injected  fluid,  such  as  the  density,  p ,  and 
viscosity,  p .  According  to  the  Buckingham  Pi  theorem,41  two  dimensionless  equations 
for  the  film  thickness  and  spreading  angle  at  the  injector  exit  can  be  obtained  in  terms  of 
those  variables, 

h/Dn  =  Apinj  /  Ao  >  Mmj  /M„,ReL,K)  ( 1 ) 

«  =  tan  Xueg/uex=f(pinj/px,pinJ/px,ReL,K)  (2) 

where  the  subscripts  inj  and  oo  denote  the  injectcd-fluid  and  ambient  conditions, 
respectively.  The  Reynolds  number  of  the  liquid  film  is  defined  as  ReL  =  pinju„L /  pjnJ. 

Among  those  parameters,  the  density  and  viscosity  ratios  arc  fixed  for  a  given  flow 
condition.  The  geometric  constant  and  injector  length  are  determined  by  the  injector 
configuration. 

A  parametric  study  is  conducted  to  investigate  the  effects  of  injector  geometry 
and  operating  condition  on  the  liquid  spray  behavior  over  a  pressure  range  of  100-200 
atm.  The  LOX  injection  and  ambient  temperatures  are  fixed  at  120  and  300  K, 
respectively.  For  reference,  the  critical  pressure  and  temperature  of  oxygen  arc  50.4  atm 
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and  154.6  K,  respectively.  Table  2  summarizes  the  injector  configurations  and  flow 
parameters  obtained  from  the  present  analysis.  Also  included  for  comparison  are  those 
predicted  by  classical  swirl-injector  theories  for  inviscid,  incompressible  flows  (denoted 
by  the  subscript  mv),1  written  in  the  following  form, 

(3) 

a,m  =  =  tan  1  pJ\-<pY  l\2-,p-2(\-:pY\  (4) 

Table  2  Effects  of  injector  geometry  and  flow  conditions  on  LOX  film  thickness  and 
spreading  angle  at  injector  exit  (Tx=  300  K,  and  Tinj  =  120  K). 


Cases 

Poo 

(atm) 

L 

(mm) 

AL* 

(mm) 

K 

m 

(kg/s) 

R  eL 
(106) 

h 

(mm) 

2  a 
(deg) 

h 

inv 

(mm) 

2a* 

(deg) 

i 

100 

25 

2.0 

3.2 

0.15 

5.4 

0.686 

73.8 

0.56 

92.2 

2 

100 

25 

2.0 

3.2 

0.10 

3.8 

0.660 

73.6 

0.56 

92.2 

3 

100 

25 

2.0 

3.2 

0.20 

7.6 

0.632 

73.2 

0.56 

92.2 

4 

100 

25 

2.0 

3.2 

0.25 

9.5 

0.611 

73.0 

0.56 

92.2 

5 

100 

25 

0.5 

3.2 

0.15 

5.7 

0.660 

76.0 

0.56 

92.2 

6 

100 

25 

4.5 

3.2 

0.15 

4.9 

0.672 

74.4 

0.56 

92.2 

7 

100 

50 

2.0 

3.2 

0.15 

11.0 

0.709 

73.4 

0.56 

92.2 

8 

100 

75 

2.0 

3.2 

0.15 

16.0 

0.738 

72.0 

0.56 

92.2 

9 

100 

100 

2.0 

3.2 

0.15 

22.0 

0.796 

69.8 

0.56 

92.2 

10 

100 

25 

2.0 

4.2 

0.15 

5.8 

0.595 

87.0 

0.50 

103.6 

11 

150 

25 

2.0 

3.2 

0.15 

5.0 

0.612 

72.4 

0.56 

92.2 

12 

200 

25 

2.0 

3.2 

0.15 

4.8 

0.586 

71.8 

0.56 

92.2 

*  distance  from  the  center  of  the  tangential  entry  to  the  injector  headend. 


where  the  coefficient  of  passage  fullness,  (p ,  is  defined  as  the  fractional  area  filled  by  the 
injected  fluid  in  the  discharge  nozzle, 


nRl  r: 


(5) 


As  will  be  further  elaborated,  both  the  liquid-film  thickness  and  spreading  angle 
are  sole  functions  of  the  injector  geometric  constant,  K  ,  and  independent  of  such 
operating  parameters  as  mass  flowrate,  in  accordance  with  classical  hydrodynamics 
theories.  The  situation  becomes  substantially  different  in  reality  due  to  the  existence  of 
viscous  loss  and  property  variations.  Since  no  distinct  interface  exists  between  the 
injected  fluid  and  gaseous  core  in  a  supercritical  flow  environment,  to  facilitate  analysis, 
the  numerically  calculated  LOX  film  boundary  is  defined  as  the  radial  position  above 
which  the  mass  flowrate  equals  the  specified  value  at  the  inlet.  This  definition  takes  into 
account  rapid  fluid  property  variations  and  provides  a  quantitative  basis  for  evaluating  the 
liquid-sheet  formation  and  breakup  processes.  Because  the  detailed  velocity  information 
can  be  easily  retrieved  from  the  present  analysis,  the  film  spreading  angle,  a  ,  is 
evaluated  based  on  its  definition,  as  the  inverse  tangent  of  the  ratio  of  the  averaged 
azimuthal  to  axial  velocity  at  the  injector  exit. 
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A  typical  simulation  employs  a  50x100  grid  for  the  injector  interior,  and  a 
150x200  grid  for  the  downstream  exterior  region.  The  grids  arc  clustered  near  the  wall 
and  the  injector  exit  to  resolve  steep  gradients  in  these  regions.  The  smallest  grid  size  in 
the  radial  direction  is  20  /urn,  which  falls  in  the  inertial  sub-range  of  the  turbulent  kinetic 
energy  spectrum  estimated  using  the  Kolmogorov-Obukhow  theory.  The  computational 
domain  is  divided  into  15  blocks,  with  each  calculated  on  a  single  processor  of  a 
distributed-memory  parallel  computer.  The  physical  time  step  is  1x10  "  ms  and  the 
maximum  CFL  number  for  the  inner-loop  pseudo-time  integration  is  0.7.  Each 
simulation  is  conducted  for  10  flow-through  times  for  the  entire  domain  (i.e.,  0.1  s)  to 
obtain  statistically  meaningful  data. 

A  grid-independence  study  was  performed  for  Case  1  as  part  of  the  validation 
procedure.  The  same  numerical  code,  injector  configuration,  and  flow  conditions  arc 
considered  using  a  fine  mesh  with  a  75x160  grid  for  the  injector  interior  and  a 
225x320  grid  for  the  downstream  exterior  region.  The  mean  grid  resolution  inside  the 
injector  is  enhanced  by  50%  in  each  spatial  direction.  Figure  3  shows  the  radial 
distributions  of  the  time-mean  temperature,  density,  and  axial  and  radial  velocity 
components  at  different  axial  locations.  Results  for  the  two  different  grid  systems  are  in 
close  agreement.  The  maximum  derivation  in  the  entire  domain  is  less  than  5%.  The 
effect  of  grid  resolution  on  the  dynamical  behavior  of  the  injector  flow  was  also 
examined.  The  maximum  relative  error  is  3%  for  both  the  length  and  speed  of  the 
hydrodynamic  instability  wave  within  the  liquid  film.  The  differences  of  the  frequency 
and  magnitude  of  the  dominant  pressure  oscillations  induced  by  flow  instabilities  are  less 
than  5%.  The  results  indicate  that  the  injector  dynamics  have  been  well  captured  and  the 
grid  system  employed  is  appropriate  in  the  present  work. 


Fig.  3  Effect  of  grid  resolution  on  radial  distributions  of  mean  temperature,  density,  and 
axial  and  radial  velocity  components  at  different  axial  locations  (Case  1:  px  - 10  MPa , 
TinJ  =  120  K ,  Tx  =  300  K ,  m  =  0. 1 5  kg / s ,  A L  =  2.0  mm ). 
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1.5.1  Injector  Flow  Dynamics 

Figure  4  shows  the  temporal  evolution  of  the  temperature  field  for  Case  1  (the  baseline 
case).  Liquid  oxygen  is  introduced  tangentially  into  the  injector.  The  swirl-induced 
centrifugal  force  prevents  the  injected  fluid  from  penetrating  into  the  center,  and 
consequently  gives  rise  to  the  formation  of  a  gaseous  core  (referred  to  as  cavity).  A  thin 
liquid  film  forms,  convects  downstream  along  the  wall  according  to  the  conservation  of 
mass  and  momentum,  and  exits  from  the  injector  as  a  nearly  conical  sheet.  Several 
different  types  of  flow  oscillations  are  observed.  The  LOX  film  is  intrinsically  unstable 
and  features  three-dimensional  hydrodynamic  instability  waves  in  both  the  longitudinal 
and  circumferential  directions. 4f’  The  circumferential  mode,  however,  can  not  be 
observed  in  the  present  axisymmetric  simulation.  The  longitudinal  instability  wave 
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Fig.  4  Temporal  evolution  of  temperature  field  (Case  1 :  px  =  10  MPa ,  7]  .  =  120  K , 

Tx  =300 K,  m  =  0.\5kg/s,  AT  =  3.2,  AL  =  2.0  mm). 

propagates  within  the  film  in  a  form  similar  to  the  shallow-water  wave  for  an 
incompressible  flow,  and  induces  oscillations  of  the  circumferential  velocity  in  both  the 
axial  and  radial  directions,  which  are  then  convected  downstream  with  the  mean  flow. 
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As  the  liquid  film  is  discharged  from  the  injector  to  the  chamber,  the  Kelvin-Helmholtz 
type  of  instability  emerges  in  the  LOX  sheet,  and  subsequently  leads  to  sheet  breakup  and 
droplet  formation.  In  addition,  the  injector  is  acoustically  connected  with  the 
downstream  chamber.  Any  acoustic  oscillation  in  the  chamber  may  propagate  upstream 
into  the  injector  and  affect  the  LOX  film  behavior  by  influencing  the  mass  flowrate  of  the 
injected  fluid. 1,6  Furthermore,  strong  acoustic  resonance  occurs  if  the  natural  frequencies 
of  the  injector  and  chamber  match  each  other.4 

Figure  5  shows  a  close-up  view  of  the  temperature  evolution  near  the  LOX  film 
for  the  baseline  case,  indicating  the  presence  of  the  longitudinal  mode  of  the 
hydrodynamic  instability.  The  wave  grows  and  develops  to  large-scale  billows  as  it 
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Fig.  5  Close-up  view  of  temporal  evolution  of  temperature  field  near  liquid-oxygen  film 
(Case  1:  =\0 MPa  TinJ=l20K  Tx  =300 K  m  =  Q.\5kgls  K  =  3.2, 

A L  =  2.0  mm ) 


convects  downstream.  The  calculated  wave  speed  is  approximately  10  m/s  in  the  region 
of  2<  x/ Rs  <9 ,  where  the  LOX  film  is  well  developed  and  its  time-mean  thickness 

varies  slowly.  A  wave  equation  characterizing  the  hydrodynamic  instability  of  a  swirling 
liquid  film  can  be  established  based  on  linearized  conservation  laws.  For  an  inviscid, 
incompressible  flow,  with  the  neglect  of  the  radial  velocity  and  the  assumption  of  an 
infinitesimal  film  thickness  compared  to  the  wave  length,  the  wave  equation  becomes6 


dt 2 


(6) 


where  £  is  the  instantaneous  displacement  of  the  film  surface,  and  rm  the  radial  distance 
between  the  injector  centerline  and  the  time -mean  film  surface,  and  uin  the  injection 
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velocity  at  the  tangential  inlet.  The  form  of  the  wave  speed  bears  a  close  resemblance  to 
that  for  shallow-water  wave  propagation,6 


ulR. 


Rl-rl 


a  =  J(^LrLK-Jr-^L) 

V  r  2  r 


(7) 


The  first  parenthesized  term  in  the  square  root  represents  the  centrifugal  force  and  the 
second  term  the  effective  thickness  of  the  liquid  film.  Equation  (7)  results  in  a  wave 
speed  of  6.6  m/s  based  on  the  values  of  um  and  rm  (11  m/s  and  1.9  mm,  respectively). 

The  classical  hydrodynamics  theory  underpredicts  the  wave  speed  due  to  the  neglect  of 
fluid  compressibility  and  viscous  effects. 

The  frequency  spectra  of  pressure  fluctuations  provide  a  more  quantitative  insight 
into  the  various  types  of  instability  waves  involved.  Figure  6  presents  the  results 
obtained  at  eight  different  positions  within  the  injector.  A  small  peak  at  14  kHz  is 
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Fig.  6  Power  spectral  densities  of  pressure  fluctuations  at  eight  different  locations  inside 

injector  (Case  1:  P.^MPa^  TL,  =120*  r„  =300 1  A  =  0.15%/S,  AT  =  3.2, 

A L  =  2.0  mm  y 
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observed  at  probes  1  -4,  resulting  from  the  recirculating  LOX  flow  between  the  tangential 
inlet  and  the  injector  head-end.  This  high-frequency  oscillation  is  confined  within  the 
upstream  region.  It  decays  rapidly  and  disappears  as  the  LOX  film  moves  downstream 
due  to  viscous  dissipation.  Two  dominant  modes  at  the  frequencies  of  0.55  and  3.15  kHz 
are  observed  at  all  the  probe  locations.  The  former  is  closely  related  to  the  longitudinal 
wave  of  hydrodynamic  instability  within  the  LOX  film.  In  the  present  case,  the  wave 
propagation  speed  is  about  10  m/s,  and  the  time  for  a  disturbance  to  travel  through  the 
LOX  film  within  the  injector  (i.e.,  from  the  tangential  entry  to  the  injector  exit)  is  on  the 
order  of  2  ms.  This  leads  to  a  characteristic  frequency  of  500  Hz  for  the  longitudinal 
hydrodynamic  instability.  The  result  shows  good  agreement  with  the  calculated 
frequency  of  550  Hz.  Figure  7  presents  the  frequency  contents  of  the  pressure  and 
velocity  fluctuations  at  Probe  6.  Both  the  axial  and  azimuthal  velocities  correlate  closely 
with  the  pressure  signal.  As  the  wave  propagates  downstream,  the  local  axial  velocity 
and  pressure  also  fluctuate  to  satisfy  the  conservation  of  the  mass  and  momentum.  The 
ensuing  variation  of  the  film  thickness  then  causes  the  azimuthal  velocity  to  oscillate  in 
accordance  with  the  conservation  of  mass  and  angular  momentum.  The  lack  of 
noticeable  oscillations  of  the  radial  velocity  in  the  low-frcquency  range  may  be  attributed 
to  the  small  thickness  of  the  LOX  film. 
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Fig.  7  Power  spectral  densities  of  pressure  and  velocity  fluctuations  at  probe  6  (Case  1  : 
px  —  10  MPa ,  Tinj  =  120  K ,  Tx  —  300  K ,  m-  0. 15  kg! s ,  K  =  3.2 ,  A L  =  2.0  mm ). 


The  harmonic  of  3.15  kHz  shown  in  the  frequency  spectral  of  Fig.  6  may  arise 
from  the  acoustic  oscillation  in  the  injector.  The  injector  configuration  considered  in  the 
present  study  can  be  treated  acoustically  as  a  quarter-wave  resonator  having  a  natural 
frequency 

f  =  c/  4(L  +  A/)  (8) 

where  c  is  the  speed  of  sound  in  the  gaseous  oxygen  core,  L  the  injector  length.  A 
correction  factor,  A/,  which  is  usually  taken  as  0.6 Rn,  is  employed  to  account  for  the 

effect  of  gaseous  oxygen  immediately  downstream  of  the  injector  exit  on  the  acoustic 
resonance  of  the  injector.  Under  the  present  flow  conditions  of  c  =  340  m/s  , 
L  —  25  mm  ,  and  Al-l.5mm,  the  resonance  frequency  becomes  3.2  kHz,  almost 
identical  to  the  observed  harmonic  in  Fig.  6.  Figure  8  shows  the  time  histories  of 
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pressure  oscillations  at  Probes  5-8.  The  low-frequency  (0.55  kHz)  hydrodynamic  wave 
decays  considerably  in  the  downstream  region,  while  the  acoustic-driven,  high-frequency 
(3.15  kHz)  instability  retains  its  magnitude  in  the  entire  LOX  film. 
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Fig.  8  Time  histories  of  pressure  oscillations  at  four  different  locations  along  liquid- 
oxygen  film  (Case  1:  -  10 MPa,  Tjnj  =  120 AT,  Tn  =300 AT,  m  =  0.15  kg/s ,  K  =  3.2, 

A L  =  2.0  mm ). 


Figure  9  shows  the  spectra  contents  of  flow  oscillations  in  the  near-field  of  the 
injector.  A  dominant  frequency  of  1.04  kHz  is  observed  at  all  probes.  This  phenomenon 
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Fig.  9  Power  spectral  densities  of  pressure  fluctuations  at  four  different  locations  in  ncar- 
field  of  injector  (Case  1:  =  1 0  MPa ,  Tinj  =  120  AT,  Tx  =300  K,  m  =  0.15  kg/ s , 

K  =  3.2 ,  A L  =  2.0  mm ). 
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can  be  attributed  to  the  precession  of  the  vortex  core  in  the  central  toroidal  recirculation 
zone  (CTRZ)  immediately  downstream  of  the  injector  exit/"  The  resultant  oscillation 
then  excites  flow  instabilities  inside  the  injector,  as  evidenced  in  the  frequency  spectral 
given  in  Fig.  6.  The  LOX  sheet,  once  it  exits  from  the  injector,  is  subject  to  the  Kelvin- 
Helmholtz  type  of  instability.  A  peak  around  3.9  kHz  (Probes  9-10)  exists  along  the 
LOX  sheet,  which  subsequently  develops  to  4.8  kHz  (Probes  11-12)  as  the  liquid  travels 
downstream.  The  vortex-shedding  frequency  can  be  roughly  estimated  using  the 
following  empirical  correlation, 4Q 


fn  =  St-U  /  9 


(9) 


where  the  Strouhal  numer,  St,  ranges  from  0.044  to  0.048  for  turbulent  flows.  In  the 
present  study,  the  mean  velocity,  U ,  is  around  10  m/s,  and  the  momentum  thickness  of 
the  conic  sheet  is  around  0.08  mm.  The  most  unstable  frequency,  fn ,  is  predicted  to  be 
5.75  kHz,  which  is  comparable  with  the  calculated  frequency. 


The  injector  flow  dynamics  are  further  explored  using  the  proper-orthogonal- 
decomposition  (POD)  technique,  which  extracts  energetic  coherent  structures  from  the 
calculated  flowfields.  For  a  given  flow  property,  f{r,t)  ,  the  POD  analysis  can 

determine  a  set  of  orthogonal  functions  <pnj  =  1,  2,  •  ••,  such  that  the  projection  of  / 
onto  the  first  n  functions 


f(r,t)  =  f(r)  +  djittyir ) 

H 


(10) 


has  the  smallest  error,  defined  as  E 


Here,  tf;(t)  represents  the  temporal 


variation  of  the  y'th  mode,  and  E()  and  |||  denote  the  time  average  and  norm  in  the  Z.2 


space,  respectively.  The  function  /  can  be  extended  to  a  vector  by  introducing  an 

appropriate  inner  product  on  F.  A  more  complete  discussion  of  this  subject  can  be 
found  in  Refs. 50  and  51. 

The  method  of  snapshots  is  implemented  to  compute  the  POD  modes.  The 
database  for  the  POD  analysis  contains  a  total  of  320  snapshots  of  the  flowfield  within 
the  injector.  The  time  interval  between  snapshots  is  50  //y,  compared  with  the  time  step 
of  5  /us  employed  in  the  numerical  simulations. 

Figure  10  presents  the  energy  distribution  of  the  POD  modes  for  the  oscillatory 
pressure  field  of  the  baseline  case.  Here  the  energy  of  the y'th  mode,  E  ,  is  defined  as 

£,.=E(||a//)^r)f)  (11) 

The  first  six  POD  modes  capture  more  than  95%  of  the  total  energy  of  the  oscillatory 
flowfield.  The  frequency  spectra  of  the  time-varying  coefficients,  a ,.(/),  of  these  modes 

are  shown  in  Fig.  11.  A  dominant  frequency  of  0.55  kHz  is  observed  for  the  first  mode, 
which  is  close  to  the  characteristic  frequency  of  the  longitudinal  hydrodynamic  instability 
wave  in  the  LOX  film.  The  second  and  third  modes  are  associated  with  the  oscillations 
resulting  from  the  precession  of  the  toroidal  recirculating  flow  downstream  of  the  injector 
exit.  The  fourth  and  fifth  modes  have  the  identical  frequency  of  3.5  kHz,  corresponding 
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to  the  natural  acoustic  oscillations  within  the  injector.  It  should  be  noted  that  the 
frequency  contents  of  the  POD  modes  deviate  slightly  from  those  presented  in  Fig.  6  (i.e., 
pressure  oscillations  in  the  LOX  film).  The  POD  analysis  is  concerned  with  the  entire 
field,  and  consequently  provides  results  in  a  volume-average  sense,  as  opposed  to  results 
at  selected  points  (Fig.  6).  The  discrepancy  becomes  more  noticeable  in  the  present 
study,  in  which  the  flowfield  is  highly  non-homogeneous  with  a  steep  density 
stratification  across  the  LOX  film  surface.  Another  factor  contributing  to  the  difference 
between  Figs.  6  and  1 1  is  the  limited  temporal  resolution  of  the  POD  analysis,  especially 
for  the  high-frequency  modes. 


Fig.  10  Energy  distribution  of  POD  modes  of  pressure  oscillations  within  injector  (Case 
1:  pm  =10 MPa,  Tinj  =\20 K ,  Tx  =300 K,  m  =  0A5kg/s,  K=  3.2,  AL  =  2.0mm). 
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Fig.  1 1  Frequency  spectra  of  time-varying  coefficients  of  first  six  POD  modes  of 
pressure  oscillations  (Case  1 :  pK  =  1 0  MPa ,  Tmj  =  120  K ,  Tx  =300  K ,  m  =  0. 1 5  kg  /  s , 
K  =  3.2 ,  AL  =  2.0  mm ). 

Figure  12  shows  the  spatial  distributions  (i.e.,  mode  shapes)  of  the  first  six  modes 
of  the  oscillatory  pressure  field  in  the  injector.  The  first  mode  shape  is  nearly  onc- 
dimensional  and  exhibits  a  decaying  distribution  with  a  maximum  at  the  head  and  a 
diminished  value  at  the  exit.  The  results  corroborate  the  intimate  relationship  between 
the  hydrodynamic  instability  wave  in  the  LOX  film  and  the  unsteady  field  in  the  gaseous 
core.  Any  hydrodynamic  disturbance  in  the  liquid  film  may  cause  variations  of  the  film 
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thickness,  which  in  turn  produces  pressure  fluctuations  in  the  free  volume  of  the  injector, 
mainly  due  to  the  volume  dilation  across  the  film  surface  induced  by  the  local 
temperature  changes  (see  Fig.  5).  The  second  and  third  modes  are  attributed  to  the 
excitation  by  the  precession  vortex  core  in  the  near-field  of  the  injector  exit.  The  fourth 
and  fifth  modes  bear  quite  similar  structures,  but  with  a  phase  difference  of  nil  in  both 
time  and  space.  This  suggests  the  existence  of  a  standing  wave  in  the  injector.  The  sixth 
mode  reveals  a  dipole-like  structure  surrounding  the  tangential  entry  near  the  injector 
headend.  The  associated  high-frequency  oscillation  is  basically  confined  in  the  upstream 
region,  and  decays  rapidly  as  the  fluctuation  propagates  downstream. 
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Fig.  12  Spatial  distributions  of  first  six  POD  modes  of  oscillatory  pressure  field  within 
injector  (Case  1:  Poo=lOMPa,  Tmj  =120  AT,  Tm  =300 K,  m  =  0.\5kg/s,  K  =  3.2, 

A L  =  2.0  mm). 

1.5.2  Mean  Flow  Properties 

Figure  13  shows  the  time-mean  fields  of  the  density,  temperature,  compressibility 
factor,  and  velocities  for  Cases  1  and  10.  The  corresponding  geometric  parameters,  K , 
are  3.2  and  4.2,  respectively.  The  critical  isothermal  surface  of  oxygen  at  Tc  =  154  K ,  as 

denoted  by  the  dashed  line,  is  also  presented,  to  provide  a  reference  for  the  LOX  film 
boundary.  The  dense  fluid  film  near  the  wall  and  gaseous  core  in  the  cavity  are  clearly 
observed.  The  LOX  sheet  mixes  promptly  with  the  warm  gas  soon  after  exiting  from  the 
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injector,  instead  of  penetrating  deeply  into  the  chamber.  Rapid  volume  dilation  and 
property  variations  take  place  when  the  LOX  mixes  with  the  surrounding  gaseous  oxygen 
and  the  local  temperature  transits  across  the  inflection  point  on  the  isobaric  line  in  the 
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Fig.  1 3  Time-mean  fields  of  density,  temperature,  compressibility  factor,  and  velocity 
components  (px  =  10  MPa ,  Tmj  =  120  K ,  Tx  =  300  K  ,  m  =  0. 15  kg/ s ,  AL  =  2.0  mm ). 


thermodynamic  p-T  diagram.  Figure  14  shows  the  radial  distributions  of  the  time-mean 
flow  properties  at  various  axial  locations  for  Case  1.  The  velocity  field  indicates  the 
presence  of  a  swirling  gaseous  flow  downstream  of  the  injector  exit,  where  the  adverse 
pressure  gradient  in  the  axial  direction  produces  a  recirculation  zone.  Owing  to  the 
relatively  low  pressure  in  this  region,  the  LOX  sheet,  which  initially  expands  outward  in 
the  radial  direction,  converges  toward  the  centerline,  until  the  radial  pressure  gradient  is 
balanced  by  the  centrifugal  force.  As  the  swirl  strength  increases  with  K  varying  from 
3.2  to  4.2,  the  film  thickness  decreases,  but  the  spreading  angle  becomes  wider.  The 
recirculation  zone  downstream  of  the  exit  is  enlarged  and  even  extended  into  the  injector. 
A  similar  situation  was  observed  by  Panda  and  Mclaughlin  2  in  their  experimental  study 
on  swirling-jet  instabilities  at  a  swirl  number  of  0.5  and  a  Reynolds  number  of  57000. 

Figure  15  shows  the  distribution  of  the  LOX  film  thickness  along  the  injector  wall 
for  Case  1,  defined  based  on  the  injected  mass  flow  rate  of  0.15  kg/s  and  local  fluid 
density  and  axial  velocity.  Three  different  flow  regimes  exist  within  the  injector.  In  the 
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Fig.  14  Radial  distributions  of  mean  velocity  components,  temperature,  density  and 
compressibility-factor  at  various  axial  locations  (Case  1:  pm=  10  MPa ,  TinJ  =120  K , 

Tx  =300 K,  m  =  Q.\5kg/s,  K  =  3.2,  AL  =  2.0mm). 


developing  region  (0<x/ Rs<2)  close  to  the  injector  head  end,  the  injected  LOX 

initially  occupies  a  substantial  fraction  of  the  free  volume,  due  to  the  small  axial  velocity. 
The  fluid  then  converges  toward  the  wall  and  forms  a  thin  film.  Downstream  of  the 
tangential  inlet  at  2<x/Rs<9,  a  stationary  regime  emerges,  in  which  the  radial 

distributions  of  the  flow  properties  vary  slowly  along  the  axial  direction.  The  film 
thickness  increases  gradually  due  to  the  decreased  axial  velocity  caused  by  viscous 
friction  and  reduced  fluid  density  resulting  from  local  heat  transfer.  Finally,  as  the  static 
pressure  induced  by  the  swirling  motion  is  converted  into  axial  momentum  close  to  the 
injector  exit,  the  axial  velocity  increases  and  the  film  thickness  decreases  in  the 
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acceleration  region  (9<x/ Rs  <10).  Also  included  in  the  figure  is  the  critical  isotherm 

of  oxygen.  The  heat  transfer  from  the  warm  ambient  gas  to  the  LOX  film  causes  the 
critical  isothermal  surface  to  regress  toward  the  injector  wall. 


Fig.  15  Liquid-oxygen  film  thickness  along  injector  wall  (Case  1 :  px  =  10  MPa , 
TlnJ=l20K,  Tw  =300 K,  m  =  0.\5kgls,  AT  =  3.2). 

Figure  16  shows  the  axial  distributions  of  the  time-averaged  axial  and  angular 
momenta  of  the  LOX  film.  The  flow  motions  in  both  the  axial  and  azimuthal  directions 
are  weak  close  to  the  headend.  The  angular  momentum  then  increases  rapidly  and 
reaches  its  maximum  at  the  trailing  edge  of  the  tangential  inlet.  The  injected  LOX  is 
driven  toward  the  wall  by  the  centrifugal  force  and  accelerates  downstream  according  to 
the  conservation  of  mass.  As  a  consequence,  the  axial  momentum  of  the  LOX  film 


Fig.  16  Axial  distributions  of  time-mean  axial  and  angular  momenta  (Case  1 : 
=10 MPa,  Tinj=mK,  T„  =300 K,  m  =  0A5kg/s,  AT  =3.2). 


increases  rapidly.  The  transport  of  azimuthal  momentum  is  more  complicated  and  can  be 
described  by  the  following  equation, 


<K 

dt 


1  dp  2  dur  Ua , 

•  = - — +  v(V2w,  +  —  — r---r)- 

pde  6  r2  dd  r1 


UUn 


(12) 
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The  first  term  on  the  right  hand  side  of  Eq.  (12)  represents  the  pressure  gradient  in  the 
azimuthal  direction,  which  vanishes  in  the  present  axisymmetrical  simulation.  The 
second  term  denotes  viscous  dissipation.  The  third  term  accounts  for  the  transfer 
between  the  radial  and  azimuthal  momenta.  This  quantity  vanishes  in  the  stationary 
region,  in  which  ur  is  practically  zero  in  the  LOX  film.  The  small  variations  of  the 
angular  momentum  downstream  of  the  tangential  entry  and  near  the  injector  exit  (the 
accelerating  region)  arise  from  the  rapid  change  of  ur  in  those  regions. 

Based  on  classical  hydrodynamics  theories,  the  axial,  «a ,  and  azimuthal,  u ^ , 

velocities  at  the  injector  exit  for  an  incompressible  inviscid  fluid  can  be  written  as 
follows,1 


yj\-2([~<p)2  !{2~(p) 

(13) 

\j2(\-<p)2 1(2- (p) 

(14) 

where  the  coefficient  of  passage  fullness,  (p,  is  a  function  of  the  geometrical  constant, 
K. 


K  =  (\-<p)yl2/<pJ<p  (15) 

Combining  of  Eqs.  (13)  and  (14)  gives  the  spreading  angle  of  the  liquid  sheet  at  the 
injector  exit,  as  defined  in  Eq.  (4).  The  LOX  film  thickness  and  spreading  angle  arc 
determined  solely  by  the  injector  geometry  and  are  independent  of  operating  conditions. 
The  situation  for  a  real  fluid,  as  treated  in  the  present  study,  however,  may  be  quite 
different,  due  to  the  many  underlying  assumptions  in  classical  theories.  Figure  17 
presents  the  radial  distributions  of  the  time-averaged  density,  axial  velocity,  and 
azimuthal  velocity  in  the  stationary  region  ( x/R  =5)  for  Case  1.  The  discrepancies 

between  the  classical  theories  and  the  results  of  the  present  work  are  obvious. 
Nonetheless,  classical  hydrodynamic  analysis  provides  much  useful  information  about 
the  flow  physics  and  can  serve  as  an  effective  guideline  for  treating  injector  flow 
dynamics. 


simulation  -  mviscid  theory 


Fig.  17  Radial  distributions  of  mean  density  and  velocity  components  in  the  stationary 
region,  x/Rs  =5  (Case  1 :  /?„  =  10  MPa,  TinJ  =  \2QK,  Tx  =300 K ,  m  =  0.15  kg/s, 

K  =3.2,  A L  -  2.0  mm). 
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1.5.3  Effects  of  Flow  Conditions  and  Geometry  on  Injector  Behavior 

The  effects  of  LOX  mass  flowrate  on  the  injector  dynamics  are  studied  in  Cases  1 
through  4  (see  Table  2).  Consistent  with  classical  hydrodynamic  theory,  the  calculated 
film  thickness  decreases  slightly  with  increasing  mass  flowrate.  The  spreading  angle, 
however,  is  almost  independent  of  the  variation  in  mass  flowrate.  Since  both  the  viscous 
and  compressibility  effects  are  neglected,  classical  theory  under-predicts  the  film 
thickness  and  over-prcdicts  the  spreading  angle. 

The  influence  of  the  tangential  inlet  position  on  the  injector  performance  is 
explored  through  Cases  1,  5,  and  6.  Figure  18  shows  the  corresponding  axial 
distributions  of  the  axial-  and  angular-momenta  of  the  LOX  film.  A  more  severe  viscous 
loss,  as  measured  by  the  relatively  smaller  increase  of  the  axial  momentum,  is  observed 
for  Case  5  ( A L  =  0.5  mm )  because  the  tangential  inlet  is  directly  attached  to  the  hcadend. 
A  similar  situation  occurs  if  the  tangential  inlet  is  placed  farther  downstream  (Case  6, 
A L  =  4.5  mm ).  Under  this  condition,  a  more  substantial  part  of  the  injected  fluid  is 
delivered  upstream  to  fill  the  volume  between  the  inlet  and  the  hcadend  and  forms  a 
recirculation  region  there.  This  is  manifested  by  the  instantaneous  streamlines  close  to 
the  headend  for  Cases  5  and  6  in  Fig.  19.  The  information  obtained  suggests  that  an 
optimal  position  for  the  tangential  inlet  exists  to  minimize  the  momentum  loss.  The  film 
thickness  and  spreading  angle  at  the  injector  exit  arc  insensitive  to  the  location  of  the 
tangential  entry  . 
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Fig.  18  Effect  of  tangential  inlet  location  on  mean  axial  and  angular  momentum 
distributions  (px  =  10  MPa ,  Tjnj  =120  K ,  Tx  -  300  K ,  m  =  0.15  kg/ s,  K  =  3.2). 
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Fig.  19  Effect  of  tangential  entry  position  on  flow  evolution  near  the  head  end 
(px  =  10 MPa,  Tmj  =  120 AT,  Tm  =300 K,  m  =  0.\5kg/s,  K  =  3.2). 

For  injectors  with  large  length-to-diamcter  (aspect)  ratios,  the  viscous  loss  along 
the  wall  exerts  a  substantial  impact  on  the  injection  process  and  henee  alters  the 
atomization  efficiency  and  spray  distribution.  Figure  20  shows  the  relative  axial-  and 
angular-momentum  losses  as  a  function  of  the  injector  length  (Cases  1,  7,  8,  and  9).  The 
losses  are  defined  as  the  relative  drops  of  the  momenta  at  the  injector  exit  in  reference  to 
the  corresponding  values  at  the  beginning  of  the  stationary  region  (i.e.,  x/Rs  =  2).  For 


Fig.  20  Mean  axial  and  angular  momentum  losses  as  functions  of  injector  length 
(pw  =\0MPa,  TlnJ=l20K,  Tm  =300K,  rh  =  0.l5  kg/s,  K  =  3.2). 

the  two  shorter  injectors  (L/ Rs  =10  and  20),  the  axial  momentum  at  the  injeetor  exit  is 

even  greater  than  that  at  the  beginning  of  the  stationary  region.  This  phenomenon  may  be 
attributed  to  the  transformation  of  the  swirl-indueed  static  pressure  to  the  axial 
momentum  within  the  acceleration  regime.  Sinee  the  azimuthal  velocity  is  influenced  by 
both  the  viscous  dissipation  and  defleetion  of  fluid  in  the  radial  direction,  the  angular 
momentum  experiences  a  more  severe  decay  than  its  axial  counterpart.  An  increase  in 
the  injector  length  narrows  the  spreading  angle.  The  simulation  results  for  Cases  1  and  7- 
9  indicate  that  an  elongation  of  the  injector  has  three  negative  impacts:  (1)  greater 
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momentum  losses;  (2)  thickened  liquid  film  and  enlarged  mean  droplet  size;  and  (3) 
reduced  sheet  spreading  angle. 

Figure  21  shows  the  LOX  film  thickness  at  the  injector  exit  as  a  function  of  the 
injector  length.  Classical  hydrodynamic  theory  under-estimates  the  film  thickness  by 
almost  30%,  especially  for  longer  injectors.  An  improved  estimation  is  given  by  the 
empirical  correlation  suggested  by  Inamura  et  al.,  in  which  viscous  loss  is  taken  into 
account,  but  the  influence  of  property  variations  under  supercritical  conditions  is 
neglected.  The  latter  accounts  for  the  approximately  1 0%  difference  between  the  present 
analysis  and  classical  theory. 


Fig.  2 1  LOX  film  thickness  as  function  of  injector  length. 

Two  different  injector  geometrical  constants  of  K=3.2  and  4.2  are  examined  in 
Cases  1  and  10,  respectively.  The  mean  flow  properties  are  presented  in  Fig.  13.  A 
greater  geometric  constant,  in  general,  results  in  a  stronger  swirling  motion,  which 
eventually  gives  rise  to  a  thinner  liquid  film  and  a  wider  spreading  angle.  Since  the 
mixing  region  downstream  of  the  injector  is  expanded,  a  better  injector  performance  is 
achieved. 

The  effects  of  ambient  pressure  on  the  swirl  injector  behavior  are  investigated  for 
three  different  pressures  of  100,  150,  and  200  atm  (i.e..  Cases  1,  11  and  12).  As  the 
ambient  pressure  increases,  the  momentum  transfer  between  the  LOX  film  and 
surrounding  gases  becomes  stronger,  and  the  momentum  loss  thus  increases. 
Consequently,  the  film  spreading  angle  decreases  from  73.8  degrees  for  100  atm  to  71.8 
degrees  for  200  atm.  On  the  other  hand,  the  elevated  pressure  tends  to  retard  the 
gasification  of  oxygen  and  reduce  the  low-density  region  within  the  LOX  film.  The  net 
result  leads  to  a  decreased  film  thickness  with  increasing  pressure.  The  overall  trend  is 
consistent  with  the  experimental  observations  of  Kim  et  al.'4 

1.6  Summary 

The  flow  dynamics  of  liquid  oxygen  in  a  pressure  swirl  injector  with  tangential 
entry  has  been  investigated  by  means  of  a  comprehensive  numerical  analysis.  The 
formulation  incorporates  real-fluid  thermodynamics  and  transport  into  the  conservation 
laws  to  render  a  self-consistent  approach  valid  for  any  fluid  thermodynamic  state. 
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The  flow  development  in  the  injector  broadly  involves  three  different  stages. 
Liquid  oxygen  is  introduced  to  the  injector  through  the  tangential  inlet,  and  occupies  a 
bulk  of  the  injector  volume  in  the  developing  region  near  the  headend.  A  thin  liquid  film 
then  forms,  due  to  centrifugal  force,  and  convects  downstream  along  the  wall  in 
accordance  with  the  conservation  of  mass  and  momentum.  The  flow  properties  and  film 
thickness  vary  slightly  in  this  stationary  region.  Finally,  the  LOX  flow  accelerates  in  the 
axial  direction  as  the  swirl-induced  pressure  converts  to  axial  momentum  close  to  the 
injector  exit. 

Several  different  types  of  instability  waves  were  identified  to  provide  direct 
insight  into  the  mechanisms  dictating  the  injector  flow  dynamics.  These  include 
hydrodynamic  instabilities  in  both  the  axial  and  azimuthal  direction  within  the  LOX  film, 
acoustic  waves  in  the  gaseous  core,  shear-layer  instabilities  in  the  LOX  sheet  downstream 
of  the  injector  exit,  and  swirling  LOX-shect  induced  flow  recirculation  near  the  injector 
exit.  The  interactions  between  those  flow  oscillations  and  their  influences  on  the  injector 
flow  behavior  were  examined  systematically  using  the  spectral  and  proper-orthogonal- 
decomposition  techniques. 

A  parametric  study  was  also  performed  to  investigate  the  effects  of  flow 
conditions  and  injector  geometry  on  the  LOX  sheet  behavior.  Results  were  compared 
with  predictions  from  classical  hydrodynamic  theories  in  terms  of  the  film  thickness, 
spreading  angle,  and  velocity  distributions.  The  rapid  property  variations  of  oxygen  fluid 
and  viscous  dissipation  play  a  decisive  role  in  determining  the  injector  characteristics. 
The  present  work  provides  a  quantitative  basis  for  optimizing  the  injector  performance. 
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Task  2 


Counterflow  Diffusion  Flames  of  General  Fluids:  Oxygen/Hydrogen  Mixtures 

Summary 

A  comprehensive  framework  has  been  established  to  study  laminar  counterflow 
diffusion  flames  for  general  fluids  over  the  entire  regime  of  thermodynamic  states.  The 
model  incorporates  a  unified  treatment  of  fundamental  thermodynamics  and  transport 
theories  into  an  existing  model,  the  DMCF  code,  for  treating  detailed  chemical  kinetic 
mechanisms  and  multi-species  transport.  The  resultant  scheme  can  thus  be  applied  to 
fluids  at  any  state.  Both  subcritical  and  supercritical  conditions  are  considered.  As  a 
specific  example,  diluted  and  undiluted  H2/O2  flames  arc  investigated  at  pressures  of  1-25 
MPa  and  oxygen  inlet  temperatures  of  100  and  300  K.  The  effects  of  pressure  p  and 

strain  rate  a  on  the  heat  release  rate  qs  ,  extinction  limit,  and  flame  structure  arc 
examined.  In  addition,  the  impact  of  cross-diffusion  terms,  such  as  the  Soret  and  Dufour 
effects,  on  the  flame  behavior  is  assessed.  Results  indicate  that  the  flame  thickness  8 , 
and  heat  release  rate  correlate  well  with  the  square  root  of  the  pressure  multiplied  by  the 
strain  rate  as  8f  ~  1  /  yfpa  and  q s  —  yfpa  .  The  extinction  limit  strain  rate  exhibits  a 
quasi-linear  dependence  with  p  .  Significant  real-fluid  effects  take  place  in  the 
transcritical  regimes,  as  evidenced  by  the  steep  property  variations  in  the  local  flowfield. 
However,  their  net  influence  on  the  flame  properties,  appears  to  be  limited  due  to  the 
ideal-gas  behavior  of  fluids  in  the  high-temperature  zone. 

2.1  Introduction 

Laminar  counterflow  diffusion  flames  provide  much  useful  information  about  the 
basic  properties  of  non-premixed  combustion.  Several  numerical  codes  incorporating 
detailed  chemical  kinetic  mechanisms  and  multi-species  transport,  such  as  the  Detailed 
Modeling  of  Counterflow  Flame  (DMCF)  code  (Darabiha  et  al.  1988  ;  Darabiha  1992), 
have  been  developed  to  study  flame  behavior  under  various  flow  conditions.  For  instance, 
the  effects  of  strain  rate  on  burning  behavior  and  flame  stability  were  examined 
systematically  for  a  variety  of  fuel/oxidizer  combinations  (Im  et  al.  1995).  Results  have 
been  applied  as  a  submodel  in  the  modeling  of  turbulent  diffusion  flames  in  the  flamelet 
regime,  in  which  the  flame  is  locally  assumed  to  bear  a  laminar  structure.  Thus,  a 
thorough  understanding  of  strained  laminar  flames  is  a  prerequisite  to  achieving 
improved  knowledge  of  more  complex  systems. 

The  majority  of  existing  studies  on  counterflow  diffusion  flames  have  been 
carried  out  at  low  and  moderate  pressures.  The  effects  of  supercritical  conditions  that 
often  occur  in  high-pressure  combustion  devices  (Yang  2000)  have  not  been  carefully 
explored.  Most  of  the  previous  studies  have  dealt  with  systems  involving  either  a 
gaseous  or  liquid-spray  fuel  against  an  air  flow.  The  influence  of  strain  rate,  inlet 
temperature,  and  radiative  heat  losses  on  flame  structures  was  numerically  investigated 
and  benchmarked  against  experimental  data.  Sung  et  al.  (1995)  have  shown  both 


31 


experimentally  and  numerically  that  the  flame  thickness  varies  inversely  with  the  square 
root  of  the  strain  rate  for  a  methanc/oxygen/nitrogen  diffusion  flame  at  one  atmosphere. 
Brown  et  al.  (1997)  studied  the  effect  of  hydrogen  dilution  by  nitrogen  for  diffusion 
flames  involving  (an  80/20  H2/N2  mixture)  and  air.  Balakrishnan  et  cd.  (1995)  examined 
the  extinction  and  ignition  limits  for  diluted  and  undiluted  H2/02  diffusion  flames  in  the 
pressure  range  of  0.25-10  bar  using  both  full  and  reduced  chemical  kinetic  schemes.  The 
critical  strain  rate  at  extinction  was  found  to  increase  rapidly  with  increasing  pressure. 
Williams  (2001)  examined  the  effects  of  transport  on  non-premixed  flame  structures  and 
extinction  characteristics,  and  observed  that  the  strain  rate  corresponding  to  the  extinction 
limit  is  sensitive  to  molecular  transport.  A  related  study  was  later  carried  out  numerically 
by  Ben  Dakhlia  et  al.  (2002)  on  diffusion  flames  involving  n-heptane/02/N2.  The  Soret 
effect  was  found  to  be  dependent  on  the  diluent  considered  (i.e.,  nitrogen  or  helium),  and 
appeared  in  the  flame  structure  and  fuel-vapor  diffusion  boundary  layer.  Juniper  et  al. 
(2003)  considered  the  counterflow  diffusion  flames  formed  by  gaseous  hydrogen 
impinging  on  a  pool  of  liquid  oxygen  and  gaseous  hydrogen  at  1  and  2  bar  using  a 
perfect-gas  law.  The  oxygen  temperature  was  set  to  90  K,  whereas  the  hydrogen 
temperature  varied  from  20  to  310  K.  Results  indicated  that  the  heat-release  rate  per  unit 
surface  area  is  proportional  to  the  square  root  of  the  pressure  multiplied  by  the  strain  rate. 
It  was  also  found  that  the  extinction  limit  strain  rate  increases  with  pressure,  a 
phenomenon  consistent  with  experimental  results  for  n-heptane/air  flames  (Niioka  et  al , 
1991). 

The  present  work  deals  with  the  effects  of  pressure  on  laminar  counterflow 
diffusion  flames.  Emphasis  is  placed  on  the  supercritical  conditions  typically  encountered 
in  high-pressure  combustion  devices  like  liquid-propellant  rocket,  diesel  and  gas-turbinc 
engines.  A  notable  example  is  provided  by  cryogenic  thrust  chambers  like  that  of  the 
Vulcain  2  engine  (Candcl  et  al.  2006)  in  which  liquid  oxygen  (LOX)  is  injected  at  a 
subcritical  temperature  of  80  K  into  a  high-pressure  environment  of  11.5  MPa.  These 
conditions  are  to  be  compared  with  the  thermodynamic  critical  temperature  and  pressure 
of  oxygen  which  are  154.8  K  and  5.04  MPa,  respectively.  Under  these  conditions,  the 
injected  LOX  heats  up  rapidly  and  its  interface  with  the  surrounding  gases  prevails  over  a 
finite  distance  from  the  injection  plane.  The  dense  core  disappears  progressively  as  mass 
is  transferred  from  this  region  to  the  jet  surroundings.  Several  experimental  (Mayer  et  al, 
1998;  Candel  et  al.  2006)  and  numerical  (Oefelein  and  Yang,  1998;  Oefclein  2006;  Zong 
and  Yang  2007)  studies  have  been  carried  out  to  characterize  the  supercritical  flame 
dynamics  of  shear  co-axial  injectors  for  hydrogen  and  methane  fuels.  Detailed  flow 
development  and  flame  stabilization  and  spreading  mechanisms  were  investigated  in  the 
near  field  of  the  injector  exit. 

Significant  real-gas  effects  featuring  steep  property  variations  take  place  when  the 
fluid  crosses  the  thermodynamic  transcritical  regime  (Yang  2000;  Zong  et  al.  2004).  In 
contrast,  the  fluid  behaves  like  an  ideal  gas  in  the  high-temperature  reaction  and  product 
zones.  Palle  et  al.  (2005)  conducted  numerical  simulations  for  unsteady  one-dimensional 
laminar  diffusion  flames  at  a  pressure  of  10  MPa  for  N2/02,  N2/Ci2H26  and  H2/02 
mixtures.  Three  different  models  of  a  one-step  reaction  in  the  form  A+r\C— >  P  were 
considered  and  cross-diffusion  terms  were  included.  The  Soret  effect  was  found  to  be 
non-negligible  for  species  with  different  molecular  weights,  especially  for  the 
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H2+(  1/2)02— >H20  reaction.  The  Dufour  effect  was  insignificant  for  all  cases.  Sohn  el  al. 
(1998)  studied  numerically  the  structures  and  acoustic  responses  of  undiluted  H2/O2 
diffusion  flames  at  pressures  up  to  10  MPa.  Both  detailed  and  reduced  chemical  kinetic 
schemes  were  employed.  The  calculated  strain  rate  at  the  extinction  limit  showed  a  quasi- 
linear  pressure  dependence.  This  result,  obtained  with  a  4-step  reaction  mechanism, 
remains  to  be  checked  with  a  more  complete  kinetic  scheme. 

The  purpose  of  the  present  study  is  to  develop  a  comprehensive  numerical  model 
capable  of  treating  counterflow  diffusion  flames  over  the  entire  range  of  fluid 
thermodynamic  states.  Both  subcritical  and  supercritical  conditions  are  examined.  Such  a 
general-fluid  approach  will  not  only  allow  us  to  explore  the  flame  behavior  under  various 
fluid  states,  but  will  also  serve  as  a  fundamental  tool  for  establishing  flame  submodels  for 
treating  turbulent  combustion  over  a  wide  range  of  pressures.  To  this  end,  a  detailed 
combustion  modeling  tool  (DMCF)  is  first  extended  by  implementing  general-fluid 
thermodynamics  theories  (Meng  and  Yang,  2003),  such  that  the  fluid  thermodynamics 
can  be  formulated  in  a  unified  manner.  A  general  balance  of  energy  is  then  derived,  and 
thermophysical  properties  are  evaluated  with  a  self-consistent  scheme  valid  for  general 
fluids  (Yang  2000).  As  a  specific  example,  this  analysis  is  employed  to  investigate  H2/O2 
and  H2/air  diffusion  flames  in  both  subcritical  and  supercritical  environments.  The 
influences  of  pressure  and  strain  rate  on  the  flame  structure  and  behavior,  as  well  as  the 
heat  release  rate,  are  examined  systematically.  In  addition,  the  Soret  and  Dufour  effects 
are  studied. 

2.2  Theoretical  Formulation 

The  physical  model  considered  herein  is  an  axisymmetric  laminar  diffusion  flame 
stabilized  near  the  stagnation  plane  of  two  opposing  streams,  as  shown  schematically  in 
Figure  1.  The  theoretical  basis  for  treating  such  a  flame  configuration  is  well  established 
for  perfect  gases  (Smooke  et  al.  1986;  Darabiha  1992;  Em  and  Giovangigli  1999). 


Flame' 


Fig.  1.  Schematic  diagram  of  a  counterflow  diffusion  flame. 


A  constant  strain  rate,  a,  defined  as  the  radial  gradient  of  the  radial  velocity,  du/dx,  at  the 
fuel  boundary,  is  assumed.  Following  the  approach  of  Meng  and  Yang  (2003),  the 
analysis  is  extended  by  incorporating  general-fluid  thermodynamics  theories  and  property 
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evaluation  schemes,  so  that  that  the  flame  behavior  over  the  entire  regime  of  fluid  states 
can  be  formulated  in  a  unified  manner. 


2.2.1  Governing  Equations 


The  balance  equations  of  mass,  momentum,  species  mass  fractions,  and  enthalpy 
for  a  counterflow  flame  with  a  fixed  strain  rate  can  be  expressed  in  the  following  forms. 
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where  p  is  the  density,  Yk,  Wk,  and  (b,  the  mass  fraction,  molecular  weight,  and  reaction 
rate  of  species  k,  respectively,  u+  =  u/ax  the  reduced  radial  velocity,  v  the  axial  velocity, 
and  p  the  viscosity.  The  partial-mass  enthalpy  of  species  k,  hk ,  is  introduced  to  account 

for  the  interactions  between  molecules  of  different  components  in  a  general-fluid  mixture 
(Meng  and  Yang  2003).  It  is  defined  by  the  total  mass  of  the  mixture,  m,  and  the  partial 
masses  of  all  constituent  components,  mk. 
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The  specific  enthalpy  of  the  mixture,  h,  thus  becomes 


(6) 


2.2.2  Thermodynamic  Properties  and  Equation  of  State  (EOS) 

Thermodynamic  properties,  such  as  internal  energy,  enthalpy,  and  constant- 
pressure  specific  heat  are  evaluated  based  on  fundamental  thermodynamics  theories. 
Each  property  can  be  conveniently  expressed  as  the  sum  of  the  ideal-gas  counterpart  at 
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the  same  temperature  and  a  departure  function  that  accounts  for  the  dense-fluid 
correction.  Thus, 
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where  the  subscript  0  refers  to  the  ideal  state  at  low  pressure.  The  departure  functions  on 
the  right-hand  sides  of  equations  (7)-(9)  arc  determined  using  an  appropriate  equation  of 
state.  In  the  present  study,  a  modified  Soave-Redlich-Kwong  (SRK)  equation  of  state 
(Soave  1972;  Graboski  and  Daubert  1978)  is  chosen  due  to  its  wide  range  of  validity  in 
modeling  the  fluid  p-V-T  behavior  and  ease  of  implementation.  It  takes  the  following 
form. 
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(10) 


where  R  is  the  universal  gas  constant,  and  W  the  molecular  weight  of  the  fluid  mixture. 
The  two  parameters,  a  and  b,  taking  into  account  the  effects  of  attractive  and  repulsive 
forces  among  molecules,  respectively,  are  calculated  with  the  following  mixing  rules. 
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where  Xk  is  the  mole  fraction  of  species  k  and  K,y  the  binary  interaction  coefficient 
(Graboski  and  Daubert  1978).  The  constants  a,  and  b,  are  determined  from  the  following 
universal  relationships, 
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where  T  and  pt  represent  the  critical  temperature  and  pressure  of  species  k , 
respectively.  The  third  parameter,  a „  is  given  as, 


a,  = 


1  +  5. 
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where  S,  is  a  function  of  the  accentric  factor,  co„ 


S, ,  =  0.48508  + 1 .55 1 7co,  -  0. 1 561 3co, 


(17) 


2.2.3  Transport  Properties 

Accurate  evaluation  of  transport  properties  is  crucial  for  the  study  of  high- 
pressure  flow  and  flame  dynamics.  To  account  for  the  continuous  variation  of  fluid 
properties  in  a  supercritical  environment  one  cannot  use  classical  techniques  which  deal 
individually  with  liquids  or  gases  .  In  the  present  study,  both  the  mixture  viscosity,  p, 
and  thermal  conductivity,  X,  are  determined  by  the  method  proposed  by  Chung  et  al. 
(1988),  established  based  on  the  Chapman-Enskog  theory  with  a  dense-fluid  correction. 
The  calculated  properties  agree  well  with  the  NIST  experimental  data  for  both  the  gas 
and  liquid  phases  (Congiunti  et  al.  2003). 

Estimation  of  the  binary  mass  diffusivity  for  a  fluid  mixture  at  high  pressures  is  a 
challenging  task,  due  to  the  lack  of  a  formal  theory  or  even  a  theoretically  based 
correlation.  Takahashi  (1974)  suggested  a  simple  scheme  for  predicting  the  binary  mass 
diffusivity  of  a  dense  fluid  by  means  of  a  corresponding-state  principle.  The  method, 
which  was  established  based  on  curve  fits  of  experimental  data  of  various  species  over  a 
wide  range  of  pressure  and  temperature,  provides  a  relatively  accurate  estimation  of  the 
diffusion  coefficient  of  a  fluid  mixture  when  the  temperature  is  greater  than  the  critical 
value.  Model  uncertainties,  however,  may  arise  when  the  reduced  temperature  of  the 
mixture  is  smaller  than  unity  because  most  of  the  data  employed  to  validate  the 
correlation  is  located  in  the  near-  or  super-critical  temperature  regime. 

2.2.4  Boundary  Conditions  and  Numerical  Method 

For  a  given  fuel/oxidizer  combination  and  flow  condition,  the  flame  structure  is 
defined  by  the  following  boundary  conditions: 
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along  with  an  additional  relation  v^o]  =  0. 

The  overall  system  of  equations  (1-4)  can  be  conveniently  combined  into  the 
following  form 
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(20) 


^+F(V)  =  0 


where  i|/  is  the  solution  vector  and  F  a  differential  operator.  At  a  steady  state,  i.e.,  F(\\i)  = 
0,  equation  (20)  and  the  associated  boundary  conditions  can  be  solved  by  means  of  a 
combination  of  time  marching  and  Newton  iteration  techniques  (Darabiha  1992).  A 
global  adaptive  grid  is  employed  to  refine  the  spatial  resolution  in  regions  with  steep 
gradients. 

2.3  Discussion  of  Results 

The  theoretical  and  numerical  framework  outlined  above  has  been  employed  to 
study  the  H2/O2  counterflow  diffusion  flames  under  a  variety  of  flow  conditions.  The 
work  consists  of  two  parts.  First,  the  flame  behavior  over  a  wide  range  of  subcritical 
pressures  is  studied  using  several  different  detailed  chemical  kinetics  schemes.  The 
impact  of  the  Soret  and  Dufour  effects  is  also  examined.  Second,  the  flame  behavior 
under  transcritical  and  supercritical  conditions  is  explored.  Emphasis  is  placed  on  the 
effects  of  pressure  and  strain  rate  on  the  flame  properties  and  heat  release  distributions. 

2.3.1  Subcritical  Pressures 

The  H2/O2  reaction  mechanism  employed  in  the  present  study  was  developed  by 
Li  et  ah  (2004).  The  scheme  is  extended  from  the  work  of  Muller  et  al.  (1998)  and 
contains  8  reacting  species  (i.e.,  H2,  O2,  H,  O,  OH,  HO2,  H2O,  and  H2O2)  and  19 
reversible  reactions.  Validated  for  a  wide  range  of  experimental  conditions  (T  e  [298- 
3000  K],  p  €  [0.25-87  bar])  for  laminar  premixed  flames  in  shock  tubes  and  flow 
reactors,  the  mechanism  has  been  implemented  with  success  in  simulating  non-premixed 
H2/air  counterflow  flames  (Fotachc  et  al.  1998  ;  Briones  and  Aggarwal  2005).  Figure  2 
shows  the  calculated  distributions  of  the  temperature  and  species  mass  fractions  at  the 
standard  condition  (i.e.,  TH  =  Ta  -  300  K,  p  =  1  bar).  The  mixture  fraction  Z  at  the 

flame  front  is  defined  (Poinsot  and  Veynante  2005)  as 
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with  5  =  8  and  Y°  =  T °  =  1 
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Two  different  software  tools  are  employed  herein. 

(i)  Designated  as  DMCF  both  thermodynamic  and  transport  properties  are 
determined  from  the  CHEMKIN  and  TRANSPORT  libraries  for  ideal  gases 
(Keeetal.  1986). 

(ii)  Designated  as  DMCF-PRF  (Partial  Real-Fluid);  thermodynamic  properties 
are  determined  based  on  the  SRK  EOS  (Meng  and  Yang  2003),  whereas 
transport  properties  are  estimated  using  the  TRANSPORT  library  (Kee  ct 
al.  1986). 

The  difference  between  the  flame  structures  as  calculated  using  the  two  approaches  is 
negligible,  as  evidenced  in  Figure  2.  The  real  fluid  thermodynamics  implemented  in  the 
present  study  is  validated  in  the  limit  of  ideal  gas.  The  flame  thickness  8/,  defined  as  the 

37 


full  width  at  half  maximum,  is  18  mm  for  a  strain  rate  of  20  s  The  temperature  reaches 
its  maximum  value  of  3050  K,  which  is  close  to  the  adiabatic  flame  temperature  of  Taj= 
3080  K,  for  the  stoichiometric  H2/O2  mixture. 


Fig.  2.  Mass  fractions  and  temperature  at  p=l  bar,  TO2=Th2=300  K  and  a=20  s'1:  Case  (i) 
is  based  on  the  ideal-gas  assumption  (Darabiha  1992)  (□□□);  case  (ii)  incorporates 
general-fluid  thermodynamics  theories  and  the  SRK  EOS  ( — ). 

Two  other  different  detailed  chemical  kinetics  schemes  were  also  tested  for 
simulating  the  flame  structure  :  the  mechanism  of  Miller  et  al.  (1982)  and  the  sub¬ 
mechanism  of  Katta  and  Roquemore  (1998).  Almost  identical  results  were  obtained, 
except  that  the  latter  predicts  a  flame  temperature  near  4400  K.  Accordingly,  only  the 
mechanism  of  Li  et  al.  (2004)  is  used  in  what  follows. 
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The  Soret  and  Dufour  effeets  are  examined.  With  the  Hirschfelder-Curtiss 
approximation  for  thermal  diffusion,  (i.e.  the  Soret  effect),  the  mass  diffusion  velocity 
can  be  expressed  in  the  following  general  form  (Kee  et  al.  1986): 
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where  0j  is  the  thermal-diffusion  ratio.  The  correction  velocity  defined  below  can  be 
evaluated  by  the  TRANSPORT  package  (Kee  et  al.  1986). 
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This  approximation  has  been  validated  by  Daguse  et  al.  (1996)  against  a  complete  kinetie 
theory  for  a  counterflow  H2/O2/N2  diffusion  flame  at  atmospheric  pressure.  Results  show 
that  the  two  approaches  are  in  good  agreement.  The  calculated  species  and  temperature 
profiles  differ  by  1%  in  location  and  the  flame  front  is  slightly  displaced.  The  present 
method  does  not  introduce  significant  error. 

Figures  3(a)  and  3(b)  show  the  temperature  distributions  of  the  H2/air  and  H2/O2 
flames  with  strain  rates  of  a  —  2000  and  8000  s  ’  ,  respectively.  The  pressure  is  set  to  1 
bar,  and  the  incoming  flow  temperature  is  300  K  for  both  fuel  and  oxidizer.  Compared  to 
the  FL/air  system,  the  larger  quantity  of  burnt  fuel  in  the  H2/O2  system  produces  a  thieker 
flame.  The  maximum  temperature  (Tmax)  is  around  3050  K  for  the  H2/O2  system  and  2250 
K  for  the  H2/air  system.  The  latter  exhibits  a  steeper  gradient  near  the  temperature 
maximum.  At  a  higher  strain  rate,  the  flame  thickness  8/-  is  reduced  and  is  found  to  vary 
with  the  inverse  square  root  of  the  strain  rate  (i.e.  with  1  lyfa  )  (Sung  et  al.  1995)  for  both 
eases.  The  maximum  temperature  Tmax  decreases  slightly  to  3030  and  1950  K  for  the 
H2/O2  and  the  H2/air  eases,  respectively.  The  impaet  of  the  Soret  effect  is  shown  by  the 
solid  symbols  in  Figure  3.  As  a  result  of  thermal  diffusion,  light  moleeules  are  driven 
towards  hot  gases  and  heavy  moleeules  move  in  the  opposite  direction,  thereby  leading  to 
a  moderate  ehange  in  the  flame  structure.  Since  hydrogen  rapidly  reacts  with  other 
species,  thermal  diffusion  plays  a  slightly  more  important  role  on  the  oxidizer  side, 
especially  near  the  flame  center  where  light  species  are  present.  The  inclusion  of  thermal 
diffusion  thus  results  in  a  slight  decrease  in  the  flame  thiekness  and  temperature  on  the 
oxidizer  side.  The  effect  becomes  more  obvious  for  H2/O2  systems  due  to  the  steeper 
temperature  gradient  in  the  flame  zone. 

Figures  4a  and  4b  show  the  temperature  profiles  for  different  pressures  at  a  fixed 
strain  rate  of  2000  s  '  for  FL/air  and  H2/O2  flames,  respectively.  For  both  systems,  the 
maximum  flame  temperature  increases  with  increasing  pressure,  but  the  flame  thickness 
exhibits  an  opposite  trend.  This  ean  be  explained  through  an  asymptotie  analysis  (Juniper 

et  al.  2003).  The  heat  release  per  unit  surfaee  area  is  proportional  to  Jpa  ,  but  the 

thermal  conductivity  remains  basically  pressure-independent.  As  a  consequenec,  the 
maximum  temperature  increases  with  pressure.  The  reduction  in  flame  thickness  with 
pressure  was  also  found  by  Law  (2006),  who  suggested  the  use  of  the  pressure-weighted 


39 


strain  rate  pa  in  correlating  the  pressure  effects.  The  flame  thickness  is  thus  proportional 
to  l/ yfpa  (instead  of  l/\fa ),  as  shown  in  the  Figure  5. 


Fig.  3.  Impact  of  the  Soret  and  Dufour  effects  on  H2/air  and  H2/O2  counterflow  diffusion 
flames  at  p=l  bar:  lines  lines  ( — ,  — )  correspond  to  diffusion  velocity  based  on  mole 
gradient,  (■, ▲)  include  the  Soret  effect,  and  (+)  represents  the  Dufour  effect. 


The  impact  of  the  Dufour  effect,  which  accounts  for  the  thermal  diffusion  induced 
by  a  concentration  gradient  ( qD ),  is  examined  by  incorporating  the  following  term  into  the 
enthalpy  equation  (4). 


c1d=-pTu 


k=l  V 


A-e;^ 


(24) 
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The  influence  on  the  flame  structure  appears  to  be  negligible,  as  evidenced  in  Figures  3b 
and  4b.  The  same  conclusion  was  reached  by  Em  and  Giovangigli  (1998)  in  their  study 
on  the  H^/air  and  CEU/air  flames. 


Fig.  4.  Impact  of  the  Soret  and  Dufour  effects  on  fVair  and  H2/O2  counterflow  diffusion 
flames  for  different  subcritical  pressures  at  a=2000  s'1:  lines  ( — ,  —  and  •  •  •  )  correspond 
to  diffusion  velocity  based  on  mole  gradient,  (■,  A,D  and  A)  include  the  Soret  effect, 
and  (+)  represents  the  Dufour  effect  on  H2/O2  flame  (p=5  bar). 

2.3.2  Supercritical  Pressures 

The  flame  behavior  in  the  supercritical-pressure  regime  of  oxygen  (p  >  5.04  MPa) 
was  investigated  by  means  of  the  general  framework  described  in  Section  2.  All  the 
thermodynamic  and  transport  properties  are  calculated  using  the  real-fluid  approach 
along  with  the  SRK  equation  of  state.  The  overall  treatment  is  designated  as  DMCF-RF 
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(full  real-fluid).  The  chemical  reaction  mechanism  is  based  on  the  detailed  scheme 
proposed  by  Li  et  al.  (2004).  Table  1  gives  the  critical  pressure  and  temperature  of  the 
species  included  in  the  reaction  scheme.  All  the  species  except  hydrogen  have  critical 
pressure  and  temperature  larger  than  those  of  oxygen,  and  thus  may  undergo 
thermodynamic  phase  transition  to  the  liquid  phase  in  regions  situated  far  from  the  flame 
zone,  depending  on  the  local  flow  conditions.  This  phenomenon,  however,  is  not  taken 
into  account  in  the  present  analysis,  because  only  small  amounts  of  these  species  are 
present  in  the  low  temperature  regions  where  fluid  condensation  may  occur. 
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Fig.  5.  Flame  thickness,  5f,  as  a  function  for  pressure  for  PL/air  and  H2/O2  counterflow 
diffusion  flames  at  two  different  strain  rates:  ( A,V)  for  a=2000  s’1  and  (0,111)  for 
a=8000  s'1. 


Table  1  Critical  pressure  Pc  and  temperature  Tc  of  species  involved  in  the  H2/02 
chemical  reaction  mechanism. 


H2  02  H  O  OH  H02  H20  H202 

Pc  (bar)  13  50.4  88.2  76  85.4  82.8  221.2  93.5 

fc(K)  33.2  154.6  404.3  367.4  443.7  487.3  647.3  544.3 


Figures  6a  and  6b  show  the  mass  fractions  of  major  and  minor  species  at  50  and 
100  bar.  The  inlet  temperature  is  fixed  at  300  K  for  both  hydrogen  and  oxygen,  and  the 
strain  rate  is  set  to  2000  s  '.  As  the  pressure  increases  from  50  to  100  bar,  the  water  vapor 
mass  fraction  remains  the  same  in  the  flame  zone,  whereas  the  flame  thickness  varies 
inversely  with  ^[p  .  Since  the  diffusive  transport  (pD)  is  pressure  insensitive,  the 
intensity  of  chemical  reactions  becomes  the  dominant  parameter  of  the  combustion 
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species  mass  fractions  toward  the  hydrogen  side.  The  situation  is,  however,  quite 
different  for  minor  species.  At  a  higher  pressure,  the  H  mass  fraction  decreases,  and  the 
relatively  inactive  HO2  species  increases  through  the  H  +  O2  +  M  —>  HO2  +  M  reaction 
pathway.  A  small  amount  of  H2O2  is  produced  through  the  HO2-H2O2  chain  mechanism. 
Figure  7  shows  the  results  for  the  oxygen  inlet  temperature  of  100  K,  with  all  the  other 
conditions  identical  to  those  in  Figure  6.  Only  slight  changes  appear  on  the  oxygen  side, 
where  a  greater  quantity  of  energy  is  needed  to  heat  up  and  then  react  with  oxygen.  As  a 
consequence,  all  species  profiles  are  shifted  toward  the  hydrogen  side.  The  real-fluid 
effects  play  a  role  through  species  transport,  inducing  a  small  increase  in  the  maximum  of 
the  intermediate  species  as  compared  with  those  shown  in  Figure  6. 


Fig.  6.  Distributions  of  mass  fractions  of  H2/02  flames  at  50  bar  (lines)  and  100  bar 
(symbols).  Th2=  To2=300  K  and  a=2000  s'1 ;  (a)  Major  species  profiles:  02  ( — • — ,  A),  H2 
(-,V),  H20  ( — ,■)  and  OH  ( — ,♦);  (b)  Radicals  profiles:  H  ( — ,D),H02  ( — • — ,0) 
and  H2O2  (-,V). 
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Fig.  7.  Distributions  of  mass  fractions  of  H2/O2  flames  at  50  bar  (lines)  and  100  bar 
(symbols).  TH2=300  K,  TO2=100  K  and  a=2000  s'1 ;  (a)  Major  species  profiles:  02  ( — * — 
,A),  H2  (— ,V),  H20  (—, ■)  and  OH  (  — (b)  Radicals  profiles:  H  (— ,D),H02  (- 
•— ,0)  and  H202  (-,V). 

The  distributions  of  the  fluid  density  and  thermophysical  properties,  including  the 
compressibility  factor  (Zc),  the  Schmidt  and  Prandtl  numbers  (Sc(k)  =  p/(p D'k)  and  Pr  = 

{\iCpyk,  respectively),  and  the  specific  heat,  were  obtained  to  explore  their  effects  on  the 
flame  structure.  Figures  8  and  9  present  the  results  for  oxygen  inlet  temperatures  of  Ta 

=  300  and  100  K,  respectively.  The  hydrogen  temperature  is  fixed  at  300  K,  and  the 
strain  rate  at  2000  s'1.  The  condition  of  50  bar  is  close  to  the  critical  pressure  of  oxygen 
and  was  chosen  to  permit  an  investigation  of  the  effects  on  flame  behavior  of  rapid 
property  variations  in  the  transcritical  regime.  In  the  case  shown  in  Figure  8,  both  the 
oxygen  and  hydrogen  temperatures  are  supercritical.  The  compressibility  factor  is  around 
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Fig.  8.  Distribution  of  thermodynamic  and  transport  properties  for  H2/02  flame.  TO2=300 
K,  Th2=300  K,  a=2000  s'1  and  p=50  bar. 

unity  over  the  entire  domain,  whereas  the  density  changes  from  67  kg/m’  on  the  oxygen 
side  to  4  kg/nr  on  the  hydrogen  side,  with  a  minimum  of  2  kg/nv  in  the  flame  zone.  The 
transport  properties  of  Pr  and  Sc(k)  vary  moderately  around  their  mean  values  in  the 
flame  zone.  In  Figure  9,  oxygen  is  injected  at  a  subcritical  temperature  of  100  K,  but 
hydrogen  remains  at  a  supercritical  state.  The  fluid  compressibility  factor  undergoes  a 
rapid  variation  from  Zc  «  0.18  to  1.0  on  the  oxygen  side  when  the  local  temperature 
increases  across  the  critical  point.  The  same  phenomenon  is  observed  for  other 
thermodynamic  and  transport  properties  due  to  the  abnormal  changes  near  the  critical 
point  of  oxygen.  In  spite  of  such  drastic  changes  of  fluid  properties  in  the  low- 
temperature  region  on  the  oxygen  side,  the  oxygen  heats  up  rapidly  and  behaves  like  a 
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perfeet  gas  before  entering  the  flame  zone.  The  flame  thickness  as  measured  by 
thetemperature  profile  thus  beeomes  almost  identical  in  Figures  8  and  9  for  7^  =100  and 

300  K,  respectively.  The  influence  of  the  oxygen  inlet  temperature  and  associated  real- 
fluid  effect  appear  to  be  limited  in  determining  the  flame  behavior.  The  transport 
parameters  (pD  or  AJCP)  are  basically  independent  of  pressure  over  most  of  the  flowfield 
and  have  values  elose  to  their  ideal-gas  counterparts. 


Fig.  9.  Distnbutions  of  thermodynamie  and  transport  properties  for  H2/O2  flame.  To2=100 
K,  Th2=300  K,  a=2000  s'1  and  p=50  bar. 

Figure  10  shows  the  temperature  distributions  for  four  different  pressures  in  the 
range  of  25-200  bar.  The  strain  rate  is  fixed  at  2000  s’1,  and  two  different  oxygen  inlet 
temperatures  of  Tn  =  300  and  100  K  are  considered.  At  a  given  strain  rate,  an  increase  in 
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pressure  leads  to  a  thinner  flame  with  a  higher  flame  temperature.  A  decrease  in  the 
oxygen  inlet  temperature  from  300  to  100  K  causes  only  a  slight  decrease  of  20  K  in  the 
flame  temperature.  The  ensuing  enlargement  of  the  flame  is  also  quite  moderate.  This 
behavior  can  be  further  explored  by  considering  the  heat  release  rate  per  unit  area,  qs . 

(25) 

t= i 


Fig.  10.  Temperature  profiles  for  H2/C>2  counterflow  diffusion  flames  at  different 
pressures.  a=  2000  s'1.  Th2=300  K  and  To2=300  K  (lines)  or  100  K  (symbols). 

Figure  11a  shows  the  heat  release  rate  as  a  function  of  the  strain  rate  for  the 
pressure  range  of  10-250  bar.  Both  the  oxygen  and  hydrogen  inlet  temperatures  are  set  to 
300  K.  For  a  given  pressure,  the  heat  release  rate  qs  first  increases  linearly  with  the  strain 
rate,  reaches  its  maximum,  and  decreases  slightly  before  extinction.  The  linear 
dependence  of  qs  on  the  strain  rate  appears  to  be  insensitive  to  pressure,  although  a  higher 

pressure  increases  the  strain  rate  and  temperature  at  the  extinction  limit.  The  overall 
result  can  be  correlated  with  the  produet  of  pressure  and  strain  rate,  pa.  Figure  lib 

shows  the  functional  relationship  of  qs  ~^[pa  . 

Figure  12  shows  the  maximum  temperature  ( Tmax )  as  a  function  of  strain  rate.  At 
low  strain  rates,  Tmax  remains  fixed.  At  high  strain  rates,  Tmax  decreases  exponentially  to 
reaeh  its  quenehing  temperature.  The  extinction  strain  rate,  aexl,  whieh  corresponds  to  the 
end  point  of  eaeh  temperature  profile  in  Figure  12,  is  approximately  proportional  to 
pressure  and  evolves  with  this  parameter  in  a  quasi-linear  manner.  Figure  13  shows  the 
flame  thickness  (8/)  in  the  pressure  range  of  50-250  bar  for  three  different  strain  rates.  As 
at  the  suberitieal  pressures  shown  in  Figure  5,  the  result  is  plotted  as  a  function  of  the  C 

parameter,  defined  as  C  =  Sfy[pa .  The  flame  thiekness  is  approximately  proportional  to 
]/  Jpa  for  the  supercritical  pressures  considered  herein. 
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Fig.  1 1.  Effeet  of  strain  rate  on  heat  release  rate  for  H2/O1  counterflow  diffusion  flames  at 
different  pressures.  Th2=300  K  and  To2=300  K. 

2.4  Conclusion 

A  comprehensive  analysis  of  laminar  eounterflow  diffusion  flames  has  been 
developed  for  general  fluids.  The  model  incorporates  fundamental  thermodynamics  and 
transport  theories,  and  ean  be  applied  to  the  entire  regime  of  fluid  thermodynamie  states. 
In  addition,  eross-diffusion  terms  sueh  as  the  Soret  and  Dufour  effects  are  included.  As  a 
specific  example,  diluted  and  undiluted  H2/O2  flames  have  been  studied  over  a  broad 
range  of  pressures  at  both  suberitical  and  supercritical  conditions.  The  effeets  of  pressure, 
strain  rate,  and  oxygen  inlet  temperature  on  the  flame  behavior  and  energy-release  rate 
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were  examined  systematically.  Results  not  only  enhance  the  fundamental  understanding 
of  the  flame  properties  under  various  flow  conditions  and  fluid  states,  but  can  also  be 
used  as  a  submodel  in  the  treatment  of  non-premixed  turbulent  combustion. 


Fig.  12.  Effect  of  strain  rate  on  maximum  temperature  for  eounterflow  diffusion  flames  at 
Th2  =300  K  and  TO2=300  K. 
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Fig.  13.  Counterflow  diffusion  flame  thickness,  8f,  as  a  function  of  pressure,  p  for  three 
different  strain  rates. 


The  calculated  flame  thickness  Sf  and  heat-release  rate  per  unit  area  were 
found  to  depend  on  the  pressure  p  and  strain  rate  a  through  the  correlations  of 
8f  ~  1  / -Jpa  and  qs  ~^Jpa  ,  respectively.  The  extinction  limit  of  the  strain  rate  evolves 
with  pressure  in  a  quasi-linear  manner.  For  H2/02  and  H2/air  mixtures,  impacts  of  the 
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Soret  and  Dufour  effects  appear  to  be  limited  with  a  pressure  increase,  due  to  enhanced 
chemical  reaction  rates  which  override  changes  in  diffusion.  Nonetheless,  the  Soret 
effect  may  still  exert  a  non-negligible  influence  on  the  extinction  strain  rate.  Significant 
real-fluid  effects  take  place  and  manifest  themselves  by  rapid  property  variations  in  the 
transcritieal  regime.  The  resultant  influence  on  the  flame  properties,  however,  is 
insignificant,  since  the  fluid  behaves  as  a  perfect  gas  when  entering  the  high-temperature 
flame  region. 
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Task  3 


Flame  Stabilization  and  Spreading  of 
Liquid  Oxygen  (LOX)  and  Methane  at  Supercritical  Conditions 

Summary 

A  comprehensive  theoretical/numerical  framework  has  been  established  to  treat 
the  turbulent  combustion  of  liquid  oxygen  (LOX)  and  methane  under  supercritical 
conditions.  The  turbulence/chemistry  interaction  is  modeled  using  the  steady  laminar 
flamelet  model.  Detailed  LOX/methane  reaction  mechanism  is  used  for  the  chemical 
reaction.  Turbulence  closure  is  achieved  by  a  large-eddy  simulation  technique.  The 
applicability  of  the  turbulent  combustion  model  has  been  carefully  accessed  by 
comparing  the  chemical  and  turbulent  time  scales  at  condition  typical  of  liquid-propellant 
rocket  engine  operation.  Results  indicate  that  the  flamelet  assumption  is  appropriate. 
The  supercritical  mixing  and  combustion  of  LOX  and  methane  in  the  vicinity  of  a  splitter 
plate  has  been  analyzed  and  the  effect  of  real-fluid  thermodynamics  on  the  cryogenic 
flame  evolution  has  been  quantified. 

3.1  Introduction 

Although  the  use  of  cryogenic  oxygen  and  hydrogen  as  propellants  in  liquid- 
propellant  rocket  engines  has  been  well  established  for  a  variety  of  launch  vehicle 
applications,  the  low  density  of  liquid  hydrogen  gives  rise  to  a  large  propellant  tank 
volume  and  booster  size.  Compared  with  hydrogen,  methane,  as  a  potential  propellant 
candidate,  possesses  many  superior  performance  and  property  characteristics.  First, 
liquid  methane  is  six  times  denser  than  liquid  hydrogen,  and  thus,  requires  a  much 
smaller/lighter  storage  vessel.  Second,  methane  offers  the  advantage  of  being  a  “soft 
cryogen”  because  its  vaporization  temperature  is  much  higher  than  that  of  hydrogen.  It  is 
easier  to  store  and  imposes  fewer  insulation  and  handling  concerns.  Third,  the  cost  of 
methane  production  is  5  to  10  times  cheaper  than  that  of  hydrogen.  Furthermore,  out  of 
common  hydrocarbon  fuels,  methane  has  the  highest  vacuum  specific  impulse  around 
370  seconds.  Consequently,  engines  using  liquid  oxygen  (LOX)  and  methane  have 
recently  attracted  considerable  interest  for  future  development  of  high-performance 
reusable  launch  vehicles  (RLV). 

A  thorough  understanding  of  propellant  injection  and  combustion  processes  is 
essential  for  a  useful  engine  design.  Very  limited  effort,  however,  was  made  to 
investigate  LOX/methane  combustion  in  liquid  propellant  rocket  engine  environments. 
Preliminary  experimental  studies  on  the  ignition  of  LOX/methane  co-axial  spray  at  a  low 
pressure  (1.5  atm)  were  carried  out  by  Cuoco  et  al.  on  the  M3  test  bench  of  DLR, 
Germany.  The  effects  of  momentum  flux  ratio  (0.1 -2.0)  and  Weber  number  (2000- 
15000)  on  the  ignition  process  were  evaluated  at  a  chamber  pressure  of  1.5  bar  and  a 
fixed  O/F  mixture  ratio  of  3.4.  The  flame  evolution  was  visualized  by  both  shadowgraph 
and  high-speed  OH  emission  images.  Two  different  ignition  scenarios,  EB1  (Explosion- 
Blowdown-reignition)  and  I  (ignition),  were  identified.  The  EB1  process  was  more  likely 
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to  occur  at  high  momentum  flux  ratios.  Compared  with  LOX/hydrogen  systems, 
LOX/methane  flame  appeared  more  easily  to  be  lifted  off  from  the  rim  of  LOX  post.  A 
phenomenon  could  be  attributed  to  the  slow  chemical  kinetics  of  methane  combustion. 

Yang  et  al.  characterized  the  stationary  atomization  and  combustion  of 
LOX/methane  and  LOX/hydrogen  sprays  at  a  fixed  chamber  pressure  of  1.5  atm  on  the 
M3  test  bench.  The  characteristics  of  central  LOX  spray  in  terms  of  intact  core  lengths 
and  droplet  numbers  were  recorded  by  high  resolution  shadowgraph  images.  The 
intensified  high-speed  OH  radical  emission  images  were  taken  to  analyze  the  flame 
anchoring  and  subsequent  spreading.  It  was  reported  that  the  effect  of  momentum  flux 
ratio  (0.22-1.67)  and  Weber  number  (1834-15000)  on  the  reactive  shear-coaxial  spray  is 
similar  to  those  observed  under  cold-flow  conditions.  The  momentum  flux  ratio  mainly 
influenced  the  primary  break-up  of  the  LOX  jet  and  thus  determined  the  length  of  the  jet 
core,  whereas  the  Weber  number  dictated  secondary  atomization  and  subsequent  flame 
spreading.  Significant  differences  exist  between  LOX/methane  and  LOX/hydrogen 
flames  at  similar  injection  conditions  defined  by  Weber  number  and  momentum  flux 
ratio.  Lifted  flames  with  greater  spreading  angle  were  only  observed  for  methane 
combustion.  Yang  et  al.  also  extended  the  analysis  to  identify  the  factors  that  dictate  the 
oscillation  of  LOX  jet  of  a  LOX/methane  spray.  The  length  of  the  jet  core  oscillated  at  a 
low-frequency  range  of  20  to  100  Hz.  The  jet  oscillation  was  prominently  influenced  by 
the  variation  of  momentum  flux  ratio.  An  increase  of  momentum  flux  ratio  gave  rise  to  a 
decrease  of  oscillation  frequency. 

Experimental  studies  on  LOX/methane  combustion  were  also  performed  at 
elevated  pressures  typically  encountered  in  operational  liquid-propellant  rocket  engines. 
Zurback  et  a/.4  '  conducted  a  preliminary  flow  visualization  of  shear  coaxial  injection  and 
combustion  of  LOX  and  methane  at  near-critical  pressures  on  the  Mascottc  test  facility  at 
ONERA,  France.  Shadowgraph  images  revealed  that  the  flame  was  attached  to  the  LOX 
post  and  spread  downstream  along  the  oxygen  jet  boundary,  wrhich  is  quite  similar  to  that 
observed  for  the  LOX/hydrogen  system.  Singla  et  al 6  performed  a  more  detailed 
experiment  of  high-pressure  oxygen  and  methane  combustion  associated  with  a  shear 
coaxial  injector.  The  temperature  of  the  oxygen  stream  remained  at  85  K,  whereas  the 
methane  stream  took  a  value  between  120  and  288  K  to  simulate  trans-  and  super-critical 
injection  conditions.  The  chamber  pressure  varied  from  4.5  to  6.0  MPa.  Emission 
images  of  excited  OH  and  CH  radicals  were  time-averaged  to  determine  the  mean  flame 
structure.  Results  indicated  that  the  flame  was  stabilized  in  the  vicinity  of  the  LOX  post 
tip  under  all  flow  conditions.  Since  vaporization  was  the  slowest  process  at  a  subcritical 
temperature,  part  of  the  unbumed  oxygen  droplets  penetrated  into  the  inner  flame.  After 
vaporization  and  mixing  with  gaseous  methane  at  the  outer  boundary  of  the  methane 
stream,  a  second  flame  with  a  greater  expansion  angle  was  formed.  Therefore,  the  flame 
featured  two  different  regions  of  light  emission,  one  surrounding  the  liquid  oxygen  jet 
and  the  other  located  close  to  the  outer  boundary  of  the  annular  methane  stream,  when 
both  LOX  and  methane  were  injected  at  a  transcritical  condition.  The  situation  changed 
if  oxygen  was  injected  at  a  subcritical  temperature  while  methane  was  gaseous.  The 
enhancement  of  turbulent  mixing  between  subcritical  oxygen  and  gaseous  methane 
streams  led  to  only  one  flame  surrounding  the  oxygen  jet.  Singla  et  al '. 7  recently 
employed  Laser-Induced  Fluorescence  (LIF)  of  OH  radical  to  identify  the  flame  structure 
of  high-pressure  shear  coaxial  LOX/methane  combustion.  High-quality  measurements 
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were  made  in  a  pressure  range  of  1.0  to  2.5  MPa.  Beyond  the  pressure  limit  of  2.5  MPa, 
additional  fluorescence  signal  originating  from  PAH  species  outside  the  flame  zone 
interferes  with  the  OH  fluorescence  and  scrambles  the  useful  information.  Over  the 
accessible  pressures,  OH  fluorescence  images  indicated  that  the  flame  develops  as  a  thin 
wrinkled  layer  spreading  near  the  liquid  oxygen  jet.  The  flame  started  from  a  location 
away  from  the  injector  lip  indicating  the  reactive  region  is  sensitive  to  the  flow 
perturbation  and  the  flame  is  less  stabilized  than  liquid  oxygen/hydrogen  flames. 

Non-intrusive  optical  measurement  techniques  were  employed  by  Santoro  and 
colleagues*  to  investigate  swirl  coaxial  injection  and  combustion  of  LOX/GCH4  at  a  near 
critical  chamber  pressures  (4.1  MPa)  and  a  mixture  ratio  of  3.0.  Detail  information  were 
obtained  on  spray  and  flame  structure  in  the  near-field.  Experimental  studies  on  high- 
pressure  LOX/methane  combustion  were  performed  by  Lux  et  al q  over  a  wide  range  of 
pressure  (4. 1-6.8  MPa)  with  the  focus  placed  on  the  identification  of  the  operating  points 
where  unstable  combustion  occurs.  Preliminary  results  indicated  that  the  propellant 
injection  velocity  ratio  and  combustion  chamber  pressure  have  substantial  effects  on  the 
near  injector  flow  and  flame  development. 

In  parallel  to  the  experimental  studies,  numerical  simulation  was  also  carried  out 
to  investigate  the  high-pressure  combustion  of  co-flowing  LOX  and  methane  through  a 
shear-coaxial  injector.  0  Several  important  features  were  identified.  The  injector 
flowfield  could  be  characterized  by  the  evolution  of  the  three  mixing  layers  originating 
from  the  trailing  edges  of  the  two  concentric  tubes  of  the  injector.  As  a  consequence  of 
the  strong  inertia  of  the  oxygen  stream  and  the  light  density  of  methane,  a  diffusion- 
dominated  flame  was  anchored  in  the  wake  of  the  LOX  post  and  propagated  downstream 
along  the  boundary  of  the  oxygen  stream.  In  addition,  an  increase  in  the  velocity  of  the 
methane  stream  significantly  enhanced  mixing  process  and  shortened  the  potential  cores 
of  both  the  LOX  and  methane  jets.  Although  the  results  are  encouraging,  many 
uncertainties  associated  with  the  numcncal  modeling  of  supercritical  combustion  still 
need  to  be  clarified.  For  instance,  turbulence/flame  interaction  has  never  been 
successfully  addressed  in  any  numerical  modeling  of  high-pressure  combustion.  A 
moment  based  approximate  reconstruction  approach  was  employed  by  Ocfelein  1  to 
account  for  sub-grid  scale  (sgs)  species  fluctuations  on  the  chemical  reaction  in  a 
simulation  of  supercritical  LOX/hydrogen  combustion.  Research  to-date,  however, 
suggested  that  the  method  is  not  reliable  when  directly  applied  to  the  filtered  chemical 
source  terms.12 

It  is  well  known  that  the  flow  Reynolds  number  increases  almost  linearly  with 
pressure.  An  increase  in  pressure  from  1  to  100  atm  gives  rise  to  an  increase  of  the 
Reynolds  number  by  two  orders  of  magnitude.  1  The  corresponding  Kolmogorov 
microscale  decreases  by  1.5  orders  of  magnitude.  One  the  other  hand,  a  recent  numerical 
study14  on  supercritical  O2/H2  counterflow  diffusion  flames  indicated  that  the  flame 
thickness  is  inversely  proportional  to  the  square  root  of  pressure  and  reduces  slower  than 
the  size  of  Kolmogorov  eddies  as  pressure  increases.  Under  supercritical  conditions,  the 
Reynolds  number  may  reach  such  a  level  that  turbulent  eddies  could  penetrate  into  the 
flame  zone  and  promote  the  occurrence  of  unsteady  phenomena  (i.e.,  local  extinction). 
This  made  the  scale  separation  assumption  inherently  employed  in  many  subgrid-scale 
(•sgs)  turbulent  combustion  models  (i.e.,  flamelet  model)  becomes  questionable.  The 
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hypothesis  of  scale  separation  implies  that  chemistry  only  occur  in  a  thin  layer  smaller 
than  the  size  of  Komogorov  eddy  while  propellants  mix  in  the  inertial  subrange  of 
turbulence. 

In  this  work,  numerical  studies  are  performed  to  investigate  supercritical  mixing 
and  combustion  of  LOX  and  methane  separated  by  a  splitter  plate.  The  theoretical  and 
numerical  frameworks  for  the  treatment  of  general-fluid  thermodynamics  have  been 
extended  to  accommodate  sgs  turbulence/flame  interactions.  Steady  laminar  flamelet 
model  with  detail  chemistry  has  been  used  as  the  turbulent  combustion  model.  The 
specific  objectives  of  the  study  are:  1)  to  establish  a  reliable  approach  to  model  turbulent 
combustion  under  supercritical  conditions;  2)  to  characterize  the  near-field  flow 
dynamics  and  combustion  processes  in  the  vieinity  of  the  splitter  plate;  and  3)  to  explore 
the  effect  of  thermodynamic  non-idealities  on  the  high-pressure  fluid  mixing  and 
dynamics  in  terms  of  turbulent  length  scales. 

3.2  Theoretical  Formulation 


3.2.1  LES  Transport  Equations  and  Sub-grid  Scale  Models 


Turbulent  closure  in  the  present  work  is  achieved  using  a  large-eddy-simulation 
technique,  in  which  large-scale  motions  are  calculated  explicitly,  whereas  eddies  with 
scales  smaller  than  the  grid  or  filter  size  are  modeled  to  represent  the  effects  of 
unresolved  motions  on  resolved  scales.  The  formulation  treats  the  Favre-filtered 
transportation  equations  of  mass,  momentum,  energy,  and  mixture  fraction  in  the 
following  conservative  form. 
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dt  dxj 
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where  overbars  and  tildes  denote  resolved-scale  and  Favre-averaged  resolved-scale 
variables,  respectively,  p,  Uj,  p,  E,  f  ,  t,j,  qi,  and  D  represent  the  density,  velocity 
components,  pressure,  specific  total  energy,  mixture  fraction,  viscous  stress  tensor,  heat 
flux,  and  diffusivity,  respectively.  The  unclosed  subgrid-scale  (sgs)  terms  in  Eqs.  (1)  - 
(4),  including  the  sgs  stresses  rsgs ,  nonlinearity  of  viscous  stress  term  Dfs ,  heat  flux 

Q.gs ,  energy  fluxes  HJgs ,  viscous  work  crsgs ,  and  scalar  flux  d>igi ,  are  defined  as, 
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A  dynamie  Smagorinsky  model  (DSM)  proposed  by  Moin  et  al.{<  is  employed  to 
close  those  sgs  terms.  The  anisotropic  part  of  the  sgs  stresses,  Eq.  (5),  is  treated  using  the 
Smagorinsky  model,  while  the  isotropic  part,  Tskfs ,  is  modeled  with  a  formulation 
proposed  by  Y oshizawa, !  6 
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The  dimensionless  quantities  CR  and  C,  are  the  compressible  Smagorinsky  model 

constants  and  are  determined  dynamically  during  ealeulation.  The  model  utilizes  the 
information  about  resolved  seales  at  the  grid-filter  level  and  at  a  coarser  test-filter  level 
with  A  >  A.  The  least-square  method  proposed  by  Lilly1  is  then  implemented  to  obtain 
the  two  model  parameters, 

=  {L9Mt)  _1  (i6) 
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The  brackets  (•)  denotes  loeal  smoothing  in  the  test  filter, 18  whieh  is  used  to  eireumvent 

the  numerical  instability  originating  from  the  dynamically  calculating  the  coefficients  of 
the  eddy-viscosity  model.  '19  A  loeal  volume-weighted  average  with  around  27  points  is 
employed.  The  explicit  forms  of  Ltj,  M0,  /?,and  or,  are, 
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where  the  hat  represents  the  test-filtered  variable.  A  Favre-filtered  variable  at  the  test- 
filter  level  is  defined  as  f  =  pf  /  p . 


The  subgrid  energy  flux  term,  H.gs ,  is  modeled  as 
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where  H  represents  the  filtered  specific  total  enthalpy.  The  turbulent  Prandtl  number, 
takes  a  conventional  value  of  0.7.  Because  mixture  fraction,  / ,  is  a  conserved  scalar,  the 

subgrid  scalar  flux,  ,  can  be  easily  modeled  by  eddy  viscosity  models  based  on 
gradient  transport  assumption  as, 


Of  =p{uj  -uj)  =  (25) 

Sc, 

A  constant  turbulent  Schmidt  number,  Sc, ,  of  0.4  is  employed.20 

It  is  noteworthy  that  accurate  modeling  of  sgs  dynamics  under  supercritical  conditions 
remains  a  challenging  task,  due  to  complications  arising  from  rapid  property  variations 
and  real-fluid  thermodynamics  and  transport.  Very  limited  effort  has  been  applied  so  far 
to  quantify  the  effects  of  real-fluid  thermodynamics  on  small-scale  turbulent  structures 
and  a  well-calibrated  sgs  model  for  supercritical  fluid  flows  is  currently  not  available. 
This  issue  will  be  addressed  in  our  subsequent  work. 


3.2.2  Subgrid-scale  Turbulent  Combustion  Model 

The  modeling  of  turbulence/chemistry  interaction  in  a  physically  meaningful 
manner  represents  a  critical  and  challenging  issue  in  the  present  study  of  supercritical 
combustion.  Figure  1,  commonly  referred  to  as  regime  diagram  for  non-premixed 

turbulent  combustion,  shows  the  variation  of  log (yfk  /u)  with  log(//  /  /),  where  /  and  it 
are  the  reference  length  and  velocity  scales  for  non-premixed  flames,  respectively,  k  and 
l,  turbulent  kinetic  energy  and  integral  length,  respectively.  1  Four  different  combustion 


58 


regimes  exist.  Laminar  flames  occur  at  low  Reynolds  numbers.  Thin  flames,  identified 
as  flamelets,  take  place  when  the  chemical  time,  tc  ,  is  smaller  than  the  Komogorov  time 

scale,  tk  .  Under  such  a  condition,  the  flame  structure  can  not  be  influenced  by  turbulent 
eddies  and  the  scale  separation  assumption  is  valid.  The  other  extreme,  thickened  flames, 
occurs  when  rc  is  greater  than  the  integral  time  scale,  r, .  In  this  regime,  turbulent  flow 

motions  greatly  enhance  the  mixing  of  fuel  and  oxidizer.  Thus,  the  combustion  process 
can  be  locally  modeled  as  well-stirred  and  is  rate  controlling.  Flamelet  structures  interact 
with  turbulence  in  the  intermediate  regime,  called  perturbed  flamelet,  when 

For  a  turbulent  diffusion  flame,  the  local  time  and  length  scales  of  flame  depend  strongly 
on  such  flow  conditions  as  the  strain  rate,  and  are  affected  by  various  unsteady  effects 
encountered.  The  flame  structures  may  vary  at  different  spatial  locations  in  the  flowfield. 


Figure  1  Regimes  for  turbulent  non-premixed  combustion  (taken  from  Lentini  ). 
Laminar  flamelet  model 

The  basis  assumption  of  the  laminar  flamelet  model  is  that  chemical  scales  arc 
shorter  than  the  Kolmogorov  scale  of  turbulent  flows.'1  Consequently,  a  turbulent  flame 
can  be  envisioned  as  a  synthesis  of  thin  reaction  zones  (i.e.,  flamelets)  embedded  in  an 
otherwise  inert  turbulent  flowfield,  and  the  inner  structure  of  the  flame  can  be  handled 
separately  from  turbulent  flow  simulations.  Instead  of  directly  treating  the  reactive  scalar 
(i.e.,  species  conservation),  the  focus  of  the  flamelet  model  is  placed  on  the  identification 
of  the  flame  surface  in  the  flowfield,  which  can  be  obtained  by  solving  the  conservation 
equation  of  the  mixture  fraction  in  a  coupled  manner  with  the  mass,  momentum,  and 
energy  equations. 

The  flame  thickness,  however,  is  typically  smaller  than  the  grid  size  employed  in 
LES  and  is  not  actually  resolved.  Therefore,  the  filtered  species  mass  fraction  of  the  ith 
species,  Y,(x,t) ,  in  each  computational  cell  should  be  evaluated  by  convoluting  the  state 
relationships,  Yff,xsl),  with  the  sgs  filtered  probability  density  function  (FDF)  of 
mixture  fraction,  P(f) ,  and  the  sgs  FDF  of  scalar  dissipation  rate,  P(xs ,)>  as  follows, 
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(26) 


?,<*.')  -  YAf.x,,)nx„W)dxA 


It  should  be  noted  that  a  statistical  independence  is  intrinsically  assumed  in  Eq.  (26) 
between  the  sgs  variations  of  mixture  fraction  and  sealar  dissipation.  The  unresolved  sgs 
fluctuation  of  mixture  fraction  is  commonly  represented  by  a  presumed  ^-shaped 
probability  density  function  (PDF)  parameterized  by  the  filtered  mixture  fraction  and  its 
sgs  variance,  which  takes  the  following  form, 


T(a)T((3) 


T(a  +  P) 


(27) 


where  f  is  the  y  -function.  The  parameters  a  and  f3  are  defined  as, 


a  -  An> :/>-.) 

(28) 

>g=(i-/)(^(17/)-i) 

/ 2 

(29) 

The  sgs  variance  of  mixture  fraction,  /"2  , 
assumption  as,2’ 

is  modeled  based  on  the  scale  similarity 

r=Khp{j-f)2ip 

(30) 

where  Kb  is  a  model  constant  chosen  as  3.  It  has  been  validated  by  many  researchers 
that  the  p  -function  pdf  provides  an  excellent  estimation  for  the  sgs  mixture  fraction 
distribution  for  non-premixed  reacting  turbulent  flows.  3’  4  For  simplicity,  the  sgs  FDF  of 
scalar  dissipation  rate,  P(xst)  »  whieh  is  typieally  assumed  to  be  a  lognormal,  is 
considered  as  a  Dirac  peak  at  the  filtered  sealar  dissipation  rate  in  the  presently  work. 
Further  investigation  on  this  issue  is  required  in  the  subsequent  studies.  The  filtered  rate 
of  sealar  dissipation,  % ,  is  modeled  based  on  an  eddy  viseosity  approach  as  suggested  by 
Girimaji  and  Zhou, 


X 


Sc  Sc,  dxj  dxj 


(31) 


The  thermochemistry  state  relation  is  established  through  a  steady-state  flamelet 
approach  featuring  a  detailed  oxygen/methane  ehemistry  with  16  speeies  and  16  chemical 
reactions. 26  This  mechanism  has  been  validated  against  experimental  data  over  a 
pressure  range  of  1  to  20  atmospheres  for  different  flame  configurations.  Calculations 
were  performed  of  counterflow  diffusion  flames  of  LOX  and  methane  at  five  different 
nominal  strain  rates  of  4000,  5000,  8000,  16000,  and  20000  s  '.  For  all  the  flame 
calculations,  the  pressure  is  fixed  at  100  atm  and  the  inlet  temperature  of  oxygen  and 
methane  is  120  and  300  K,  respectively.  It  has  been  noted  by  Ribert  et  al.'4  that 
eryogenie  flame  is  extremely  resistant  to  flow  strain.  A  distinction  stain  rate  of  i.76xio\ 
which  is  at  least  one  order  greater  than  the  maximum  strain  rate  expeeted  in  the  jet  flame 
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within  a  typical  liquid-propellant  rocket  engine  is  obtained  under  the  present  flow 
configuration.  Therefore,  no  flame  solution  at  extinction  state  is  considered  for  the 
establishment  of  flame  library.  For  each  strain  rate,  the  scalar  dissipation  rate  is  a 
function  of  the  mixture  fraction.  To  be  consistent  with  the  flamelet  assumption,  the 
corresponding  scalar  dissipation  rate,  x  >  of  each  flame  solution  is  evaluated  at  the 
stoichiometric  location.  The  effect  of  flow  strain  on  the  flamelet  solution  is  indictated  in 
Figure  2.  The  increase  of  strain  rate  enhances  the  heat  and  species  transport  close  to  the 
flame  front,  which  leads  to  the  deceases  of  inner  layer  temperature  and  flame  thickness. 
The  production  of  OH  radial  is  also  retarded.  The  effect  of  strain  rate  on  the  thickness  of 
a  counterflow  diffusion  flame  of  gaseous  methane  and  liquid  oxygen  under  supercritical 
pressure  was  systematically  investigated  in  Ref.  14.  At  a  constant  ambient  pressure,  the 
flame  thickness  is  inversely  proportional  to  the  square  root  of  strain  rate.  To  facilitate  the 
establishment  of  the  flamelet  library,  each  flamelet  solution  obtained  has  been  expressed 
as  a  function  of  mixture  fraction,  as  shown  in  Figure  3.  The  solutions  are  then  integrated 

based  on  Eq.  (26)  and  tabulated  as  functions  of  f ,  and  f"2  .  The  present  table  has 

five  support  points  in  x  direction  and  is  resolved  by  51x51  nodes  in  fx  f2  space. 


Figure  2  Effect  of  strain  rate  on  flame  thickness  of  a  liquid  oxygcn/methanc  laminar 
diffusion  flame. 


Figure  3  Effect  of  strain  rate  on  the  behavior  of  laminar  diffusion  flames. 
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The  calculated  filtered  mixture  fraction,  mixture  fraction  variance,  and  scalar  dissipation 
rate  from  LES  simulation  are  employed  to  determine  the  appropriate  entry  of  the  flamelet 
table. 


3.2.3  Real  Fluid  Thermodynamics  and  Transport 


An  essential  prerequisite  of  any  realistic  analysis  of  supercritical  fluid  dynamics 
lies  in  the  establishment  of  a  unified  property  evaluation  scheme  capable  of  treating 
thermodynamic  properties  of  pure  substances  and  mixtures  over  the  entire  fluid 
thermodynamic  states  of  concern  from  compressed  liquid  to  dilute  gas.>  At  elevated 
pressures,  thermodynamic  models  normally  employed  to  represent  ideal-gas  behavior 
may  encounter  many  difficulties.  From  the  microscopic  point  of  view,  the  intermolecular 
mean  free  paths  reduce  with  increase  of  pressure,  and  consequently  the  molecular  volume 
and  intermolecular  forces  are  no  longer  negligible.  For  convenience,  each  property  can 
be  expressed  as  the  sum  of  an  ideal-gas  property  at  the  same  temperature  and  a 
thermodynamic  departure  function,  which  takes  into  account  the  dense-fluid  correction. 
For  example,  the  specific  internal  energy  of  a  general  fluid  mixture  can  be  expressed  as 


e(T,p)  =  e0(T)  + 
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(32) 


The  subscript  0  refers  to  a  reference  thermodynamic  state  in  the  limit  of  an  ideal  gas. 

Because  the  intermolecular  interactions  between  molecules  of  different 
components  in  a  mixture  may  differ  significantly  from  those  between  the  same  types  of 
molecules  in  a  pure  substance,  the  concepts  of  partial  properties  should  be  implemented 
to  estimate  the  properties  of  each  component  under  high  pressures.  To  facilitate 
numerical  implementation,  the  concepts  of  partial-mass  and  partial-density  properties 
were  introduced  by  Meng  and  Yang.  y  The  partial-density  internal  energy  is  defined  as, 


e=(^)r 


(33) 


where  pt  is  the  density  of  the  ith  species  of  a  fluid  mixture.  The  indices,  /  and  j ,  range 
from  1  to  N. 

An  equation  of  state  valid  for  a  broad  range  of  fluid  states  is  required  to  represent 
the  fluid  P-V-T  behavior  of  a  supercritical  fluid  mixture.  A  modified  Soave-Redlick- 
Kwong  (SRK)  equation  of  state  is  adopted  in  the  present  work  because  of  its  wide  range 
of  validity  and  ease  of  implementation.  It  takes  the  following  form, 

_pRJ__aa_P_ 

W-bp  W  ( W  +  bp ) 

where  Ru  is  the  universal  gas  constant  and  W  the  molecular  weight.  The  parameters  a 
and  h  account  for  the  effects  of  attractive  and  repulsive  forces  between  molecules, 
respectively.  The  third  parameter  a  includes  the  critical  compressibility  factor  and 
acentric  (size-shape)  interactions  between  molecules. 
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Based  on  the  SRK  equation  of  state,  the  specific  internal  energy  of  a  general  fluid 
mixture  is  derived  as, 


T2  .daa/T.  , 

e(T,p)  =  ln( 


W  +  bp 
W  +  bpJ 


(35) 


where  Yj  is  the  mass  fraction  of  the  ith  constituent  species  of  a  mixture.  Substitution  of 
Eq.  (35)  into  Eq.  (33)  yields  the  partial-density  internal  energy  of  species  i. 
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where  N  is  the  number  of  species  in  the  mixture  and  xt  is  the  mole  fraction  of  species  i. 
The  partial-density  enthalpy,  ,  can  be  easily  written  as, 


*-*+(*S- 

The  constant-volume  heat  capacity  can  be  obtained  as, 


bW  dT 


W  +  bp0 


(37) 


(38) 


Details  of  the  derivation  and  implementation  of  real-fluid  thermodynamic  properties  can 
be  found  in  Ref.  29. 


Accurate  evaluation  of  transport  properties  is  also  crucial  for  the  modeling  of 
high-pressure  flow  and  flame  dynamics.  To  aceount  for  the  continuous  variation  of  fluid 
properties  in  a  supercritical  environment,  one  cannot  use  classical  techniques  which  deal 
individually  with  liquids  or  gases.  In  the  present  study,  both  the  mixture  viscosity,  p,  and 
thermal  conductivity,  X,  are  determined  by  the  method  proposed  by  Chung  et  al„  30 
established  based  on  the  Chapman-Enskog  theory  with  a  dense-fluid  correction.  The 
calculated  properties  agree  well  with  the  NIST  experimental  data  for  both  the  gas  and 
liquid  phases.3' 


3.3  Numerical  Method 


The  theoretical  formulation  outlined  above  requires  a  robust  computational 
scheme,  due  to  the  numerical  stiffness  arising  from  rapid  flow  property  variations  and 
wide  disparities  of  the  characteristic  time  and  length  scales  involved.  To  this  end,  a 
unified  treatment  of  general  fluid  thermodynamics,  based  on  the  concepts  of  partial-mass 
and  partial-density  properties  is  established  and  incorporated  into  a  preconditioning 
scheme,29  3  .  All  the  numerical  properties,  including  the  preconditioning  matrix,  Jacobian 
matrices,  and  eigenvalues,  are  derived  directly  from  fundamental  thermodynamics 
theories,  rendering  a  self-consistent  and  robust  algorithm.  The  numerical  formulation  can 
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accommodate  any  equation  of  state,  and  is  valid  for  fluid  flows  at  all  speeds  and  at  all 
fluid  thermodynamic  states.  Further  efficiency  is  achieved  by  employing  temperature 
instead  of  enthalpy  as  the  primary  dependent  variable  in  the  preconditioned  energy 
equation.  The  procedure  eliminates  laborious  iterations  in  solving  the  equation  of  state 
for  temperature,  and  consequently  facilitates  the  load  balance  among  computational 
blocks  in  a  distributed  computing  environment.  The  resultant  scheme  is  highly  efficient 
and  suitable  for  parallel  processing. 

The  numerical  framework  employs  a  density-based,  finite-volume  methodology 
along  with  a  dual-time-step  integration  technique.  Temporal  discretization  is  obtained 
using  a  second-order  backward  differencing  scheme,  and  the  inner-loop  pseudo-time  term 
is  integrated  with  a  four-step  Runge-Kutta  scheme.  Spatial  discretization  is  achieved 
with  a  fourth-order,  central-difference  scheme  in  generalized  coordinates.  A  nine-point 
stencil  is  employed  to  evaluate  the  convective  flux  in  each  spatial  direction  to  improve 
the  spectral  resolution  of  small-scale  turbulence  structures.  Fourth-order  scalar 
dissipation  with  the  coefficient  =0.001  is  applied  to  ensure  numerical  stability  with 
minimum  contamination  of  the  solution. 

The  overall  accuracy  of  the  present  scheme  within  the  context  of  LES  was 
carefully  assessed  based  on  the  decay  of  the  kinetic  energy  of  isotropic  turbulence. 
Calculations  were  performed  for  an  isotropic  turbulent  flow  in  a  cubic  box  of  a  non- 
dimensional  width  of  2n.  The  experimental  results  of  Comet-Bellot  and  Corrsin  (CBC) 
were  selected  as  the  benchmark  with  an  initial  Taylor  Reynolds  number  of  80  on  a 
32  x  32  x  32  grid1  Figure  4  shows  the  temporal  evolution  of  resolved  turbulent  kinetic 
energy.  The  results  with  a  subgrid-scalc  model  agree  well  with  the  experimental  data, 
whereas  a  slower  decay  of  TKE  is  observed  when  the  sgs  model  turned  off.  This 
indicated  that  the  present  numerical  algorithm  offers  a  reasonable  predictive  capability  in 
terms  of  its  numerical  dissipation. 


Figure  4  Decay  of  the  resolved  turbulent  kinetic  energy. 

In  the  present  study  of  LOX/mcthane  combustion,  exceedingly  large  property 
gradients  exist  in  the  near  field  of  the  injector  where  the  fluid  density  varies  by  a  factor  of 
1 000  over  a  fraction  of  the  splitter  plate  (0.3  mm).  A  second-order  scalar  dissipation  with 
a  total-variation-diminishing  switch  developed  by  Swanson  and  Turkel34  was  thus  applied 
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to  suppress  numerical  oscillations  in  regions  with  steep  property  variations.  A  small 
dissipation  coefficient  of  0. 1 25  was  cautiously  selected  to  minimize  the  contamination  on 
the  numerical  solution. 

Finally,  a  multiblock  domain  decomposition  technique,  along  with  static  load 
balance,  is  employed  to  facilitate  the  implementation  of  parallel  computation  with 
message  passing  interfaces  (MPI)  at  the  domain  boundaries.  The  parallelization 
methodology  is  robust  and  the  speedup  is  almost  linear. 

3.4  Computational  Domain  and  Boundary  Conditions 

Figure  5  shows  the  physical  model  under  consideration.  Co-flowing  methane 
(upper)  and  oxygen  (lower)  are  injected  at  the  inlet  and  separated  by  a  0.3  mm  splitter 
plate.  The  thickness  of  the  splitter  plate  is  chosen  to  be  comparable  with  the  LOX  post 
thickness  of  the  shear  coaxial  injector  typically  employed  in  liquid-propellant  rocket 


Figure  5  Schematic  diagram  of  the  computational  domain  employed  for  the  analysis. 

engines.  The  computational  domain  includes  the  inlet  and  a  downstream  region  that 
measures  \26  and  48/>  in  the  axial  and  vertical  directions,  respectively.  The  velocity 
profile  for  a  fully  developed  turbulent  pipe  flow  is  assumed  at  the  inlet  for  both  streams. 
Turbulence  is  provided  by  superimposing  broad-band  white  noise  onto  the  mean  velocity 
profile.  The  disturbances  are  generated  by  a  Gaussian  random-number  generator  with  an 
intensity  of  5%  of  the  mean  quantity.  At  the  downstream  boundary,  extrapolation  of 
primitive  variables  from  the  interior  may  cause  undesired  reflection  of  waves  propagating 
into  the  computational  domain.  Thus,  the  non-reflecting  boundary  conditions  based  on 
the  method  of  characteristics  (MOC)  are  applied,  along  with  the  specification  of  a 
reference  pressure.  At  the  lower  and  upper  boundaries,  the  pressure  and  temperature  are 
specified  as  the  ambient  values.  The  velocity  components  are  extrapolated  from  the 
interior.  Finally,  the  no-slip  adiabatic  condition  is  enforced  along  the  solid  walls. 

3.5  Results  and  Discussion 

The  laminar  flamelet  model  is  employed  to  study  supercritical  mixing  and 
combustion  of  LOX  and  methane  separated  by  a  splitter  plate.  Liquid  oxygen  and 
gaseous  methane  are  injected  at  temperatures  of  122  and  300  K,  respectively.  The  bulk 
velocities  of  the  two  streams  are  1 0  and  60  m/s,  respectively.  The  momentum  flux  ratio 
of  the  two  streams  is  2.7,  which  is  typical  of  the  shear  coaxial  injector  for  liquid- 
propellant  rocket  engines.  The  chamber  pressure  of  10  MPa  is  well  above  the  critical 
pressure  of  oxygen  (5.04  MPa)  and  methane  (4.60  MPa).  The  critical  temperatures  of 
oxygen  and  methane  are  154  and  191  K,  respectively.  Therefore,  the  methane  stream 
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remains  at  a  supereritieal  state,  while  the  LOX  stream  is  injected  at  a  transeritieal 
temperature.  The  flow  Reynolds  number  evaluated  based  on  the  properties  of  methane 
stream  and  the  thiekness  of  the  splitter  plate  is  7.6  x  104 . 

The  computational  grid  consists  360x150  points  along  the  axial  and  vertical 
directions,  respectively.  The  grids  are  clustered  in  the  shear  layers  and  near  the  splitter 
plate  to  resolve  rapid  property  variations  in  those  regions.  The  smallest  grid  size  in  the 
radial  direction  is  6  pm ,  whieh  well  falls  in  the  inertial  sub-range  of  the  turbulent  kinetie 
energy  spectrum  estimated  using  the  Kolmogorov-Obukhow  theory.  The  computational 
domain  is  divided  into  44  bloeks,  with  eaeh  ealeulated  on  a  single  processor  of  a 
distributed  computing  facility.  The  physical  time  step  for  all  calculations  is  100  nano¬ 
seconds  and  the  maximum  CFL  number  for  the  inner-loop  pseudo-time  integration  is 
fixed  at  0.7.  For  all  of  the  models  considered,  the  simulations  are  first  eondueted  by  0.4 
ms  to  allow  the  initial  transients  to  pass  through  the  computational  domain.  The 
instantaneous  flow  properties  are  then  time  averaged  over  an  additional  1.2  ms,  around  5 
flow-through  times,  to  gather  the  mean  properties. 

To  explore  the  applicability  of  different  turbulent  combustion  models  under  the 
present  simulation  conditions,  a  comparison  between  the  ehemical  time  seale  and  the 
charaeteristie  turbulent  time  scale  is  made  first.  The  ehemical  time  is  estimated  based  on 
the  method  suggested  by  Peters3^  as. 


t  = 


zlQ-zp 


(40) 


Zsl  is  the  value  of  mixture  fraction  at  the  stoiehiometrie  condition  and  is  0.2  for 
oxygen/methane  combustion.  The  extinction  sealar  dissipation  rate,  %q ,  is  ealeulated 
from  a  eounterflow  diffusion  flame  at  a  near  extinction  state,  which  takes  the  value 
1 .56  x  105  s'1 .  The  Kolmogorov  time  scales  are  given  by, 
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where  the  integral  length  seale  is  approximated  by  the  thiekness  of  the  splitter  plate.  The 
turbulent  kinetic  energy,  k ,  viscous  dissipation  rate,  £ ,  and  kinetic  viscosity,  v ,  are 
evaluated  from  the  simulation  results.  For  instance,  viscous  dissipation,  £ ,  is  calculated 

based  on  its  definition  as,  £  =  v[Vu’  +  Vw'r] :  Vn' .  Figure  6  shows  the  ratios  of  ehemieal 
to  the  integral  and  Kolmogorov  time  scales  at  different  axial  locations.  The  turbulent 
time  seales  are  estimate  based  on  the  simulation  of  laminar  flamelet  model. 

The  weak  turbulent  kinetic  energy  in  the  developing  region  elose  to  the  splitter 
plate  renders  a  Kolmogorov  time  seale  that  is  considerably  greater  than  the  chemical 
ones.  Although  the  two  turbulent  seales  experience  significant  reductions  through  the 
transition  regime,  the  ratio  of  the  Kolmogorov  to  the  ehemical  time  remains  greater  than 
unity  throughout  the  entire  flowfield.  The  result  confirms  it  is  appropriate  to  use  flamelet 
model  under  the  present  simulation  conditions.  Cautious,  however,  still  need  to  be 
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exerted  owing  the  uncertainties  associated  with  numerical  model  and  the  intrinsic 
approximation  in  models  employed  to  evaluate  the  chemical  and  turbulence  time  scales. 


Figure  6  Axial  distributions  of  the  ratio  of  chemical  time  to  the  integral  and  Kolmogorov 
time  scales. 

Figure  7  presents  the  time-averaged  fields  of  temperature  and  mixture  fraction. 
Two  high-temperature  regions  are  observed  in  the  mean  temperature  field.  A  small 
intensive  reaction  zone  appears  immediately  behind  the  splitter  plate  and  the  main  flame 
starts  around  2  mm  downstream.  The  existence  of  a  relative  low-temperature  region 
between  those  two  reaction  zones  may  be  attributed  to  the  development  and  evolution  of 
large-scale  vortices  in  the  near-field  and  the  associated  high  flow  strain.  The 
phenomenon  also  indicates  that  the  effect  of  scalar  dissipation  on  the  finite  rate  chemistry 
has  been  well  accommodated  in  the  laminar  flamelet  model. 


x/S  xJ§ 


Figure  7  Time  averaged  temperature  and  mixture  fraction  fields 

Figure  8  shows  the  vertical  distributions  of  the  time-averaged  temperature, 
density,  and  axial  velocity  at  different  axial  locations.  Severe  density  variations  occur 
across  the  LOX  stream  boundary,  although  the  distinction  between  liquid  and  gas  phases 
diminishes  under  supercritical  pressures.  In  axial  velocity  field,  negative  mean  velocities 
appear  in  the  region  immediately  behind  the  splitter  plate  indicating  the  existence  of 
recirculation  eddies.  The  strong  recirculating  backflow  in  the  vicinity  of  the  splitter  plate 
may  act  as  a  hot  product  tank  providing  the  energy  to  ignite  incoming  propellants. 

Figure  9  shows  instantaneous  field  of  velocity  vector  and  mass  fraction  of  OH  and 
CO  in  the  vicinity  of  the  splitter  plate.  The  red  contour  in  the  velocity  vector  fields 
corresponds  to  the  isocontour  line  of  Zst=0.2,  which  represents  the  flame  boundary.  The 
large-scale  vortices  emerging  from  the  upper  comer  of  the  splitter  plate  facilitate  the 
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mixing  between  incoming  methane  stream  and  the  hot  products.  Driven  by  those  strong 
vortices,  a  relative  weak  recirculation  flow  region  forms  close  to  the  lower  comer  of  the 
splitter  plate  and  carries  the  oxygen  rich  product  toward  the  methane  stream.  Unlike 


Figure  8  Distributions  of  mean  temperature,  density,  and  axial  velocity  in  vertical 
direction  at  various  axial  locations 


Figure  9  Snapshots  of  velocity  vector,  and  mass  fractions  of  CO  and  OH  predicted  by 
laminar  flamelet  model. 
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LOX/hydrogen  combustion,  in  which  the  flame  is  anchored  very  close  to  the  LOX  stream 
because  hydrogen  is  highly  diffusive,  the  ignition  occurs  when  oxygen  and  methane  rich 
hot  products  meet  at  some  locations  between  two  propellant  streams.  Consequently,  the 
intensive  reaction  zone  starts  from  a  location  very  close  to  the  upper  comer  of  the  splitter 
plate,  as  indicated  in  the  OH  radical  or  CO  mass  fraction  field.  Since  the  flame  is 
anchored  at  a  location  very  close  to  the  high-speed  methane  stream,  the  flame  may  be 
sensitive  to  flow  perturbations  and  could  easily  be  lifted  off  from  the  splitter  plate.  The 
stability  of  liquid  oxygen  methane  flames  is  weaker  than  liquid  oxygen/hydrogen  flames. 


Figure  10  shows  the  instantaneous  fields  of  temperature,  density,  mixture  fraction, 
vorticity,  and  mass  fractions  of  O2,  CH4,  CO2,  and  H2O  in  the  near  field.  A  diffusion- 
dominated  flame  emanates  immediately  from  the  splitter  plate  and  propagates 
downstream  along  the  surface  of  the  LOX  stream.  Large  scale  vortices  emerge  from  the 
upper  comer  of  the  splitter  plate  and  greatly  enhance  the  mixing  of  methane  and  hot 
products.  Intense  property  variations  occur  as  the  fluid  state  transfer  from  cryogen 
oxygen  to  hot  product.  For  example,  the  fluid  density  varies  from  10  to  1000  kg/m  '  over 
a  very  small  interval  close  to  the  splitter  thickness  (0.3  mm).  Besides  transition  of  fluid 
states,  the  density  stratification  is  further  enhanced  by  anomalies  variations  of  specific 
heat  and  thermo  conductivity  across  the  LOX  stream  boundary  under  such  a  supercritical 
condition.  0  Consequently,  the  exceedingly  large  density  gradients  in  the  region 
surrounding  the  LOX  stream  approach  the  behavior  of  a  contact  discontinuity."  Those 
regions,  which  behave  similarly  to  a  rigid  flat  plate  within  the  flowfield  in  horizontal 
direction,  increase  the  turbulence  anisotropy  at  large  scales. 1f>  The  integral  vortices 
emerging  from  the  upper  comer  of  the  splitter  plate  evolve  in  a  manner  analogous  to  that 
produced  by  a  backward-facing  step  and  mainly  reside  on  the  lighter  methane  stream. 
The  denser  oxygen  stream  is  less  influenced.  The  growth  of  mixing  layer  between  the 
two  propellant  streams  is  retarded.  Owing  to  the  substantial  damping  effects  on  the 
turbulent  flow  in  the  vertical  direction,  the  development  of  flame  was  also  restricted. 

To  further  investigate  the  effect  of  density  stratification  on  the  cryogenic  flame 
dynamics,  the  spatial  correlations  of  mixture  fraction  have  been  calculated  at  different 
stoichiometric  mixture  locations.  The  spatial  correlation  coefficient  of  mixture  fraction  at 
spatial  location  xj  takes  the  form, 


R(x,)  = 


/'(*,- Ax)  x /'(*,.  + Ax) 
f'(x;)2 


(43) 


where  the  over  bar  denotes  a  time-averaging  procedure  and  Ax  the  spatial  increment  in 
either  axial  or  vertical  direction.  In  the  present  study,  six  pairs  of  axial  and  vertical 
mixture  fraction  signals  with  separation  distances  from  0.3£  to  1 .8<?  in  a  0.3<?  increment 
were  obtained  at  each  spatial  location  of  interest.  Figure  1 1  presents  the  spatial 
correlation  coefficients  of  mixture  fraction  in  both  the  axial  and  vertical  directions  for 
x/S  =  12,  32,  50.  The  spatial  correlation  distributions  in  both  directions  become  wider 
as  the  sampling  location  moves  downstream.  This  implies  that  the  integral  length  scale 
grows  up  with  the  development  of  large  scale  vortices.  Compared  with  the  correlation  in 
the  axial  direction,  the  correlation  in  the  vertical  direction  exhibits  much  narrower 
distributions  at  all  three  locations.  The  large  density  gradient  region  leads  to  strong 
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anisotropy  of  the  turbulence  within  the  flowfield.  It  behaves  like  a  solid  wall  that 
amplifies  the  horizontal  turbulent  fluctuation  but  damps  the  vertical  ones.  Since  eddies  of 
integral  length  scales  are  squashed  in  the  vertical  direction,  the  length  scale  in  vertical 
direction  is  much  shorter  than  horizontal  one. 
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Figure  10  Snapshots  of  distributions  of  temperature,  density,  mixture  fraction, 
vorticity,  and  mass  fractions  of  O2,  CH4,  CO2,  CO,  H2O,  and  OH  radical. 


Figure  11  Spatial  correlation  of  mixture  fraction  in  both  axial  and  vertical  directions  at 
three  different  axial  locations  along  the  mean  stoichiometric  mixture  fraction  contour. 
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Figure  12  shows  the  power  spectral  densities  of  pressure  fluctuations  at  six 
different  locations.  A  dominant  frequency  of  1 1 .2  kHz,  corresponding  to  the  vortex 
shedding  frequency  downstream  of  the  splitter  platter,  is  obtained.  Owing  to  vortex 
interactions,  higher  harmonics  exist  at  upstream  locations  (i.e.,  probes  1-4).  The 
influence  of  high  frequency  oscillations,  however,  diminishes  rapidly  as  vortices 
propagate  downstream  (i.e.,  probes  5  and  6).  It  is  noteworthy  that  the  amplitude  of  the 
dominant  pressure  oscillation  at  the  probes  surrounding  the  methane  stream  is  higher  than 
those  obtain  around  LOX  jet  boundary.  This  further  collaborates  that  the  formation  of 
large  density  gradient  region  close  to  the  LOX  stream  inhibits  the  evolution  of  large  scale 
vortices. 


r  (K.) 

II 20 

995 

1870 

2795 

3620 


a 

"Cu 


Frequency,  kHz 


Frequency,  kHz 


Frequency,  kHz 


Figure  12  Frequency  spectra  of  pressure  oscillations  at  six  different  probe  locations. 

Summary 

The  unified  theoretical/numerical  frameworks  for  the  treatment  of  general  fluid 
thermodynamics  have  been  extended  to  accommodate  turbulence/flame  interactions.  A 
laminar  flamelet  model  has  been  implemented,  which  could  accommodate  the  detail 
LOX/methane  combustion  mechanism  and  the  effect  of  flow  strain  on  the  finite-rate 
chemistry.  The  applicability  of  the  turbulent  combustion  model  has  been  accessed 
through  a  careful  comparison  of  chemical  to  characteristic  turbulent  time  scales  under  a 
condition  typical  of  liquid-propellant  rocket  engine  operation.  The  results  indicate  that 
the  flamelet  assumption  is  appropriate  and  supercritical  combustion  process  is  mixing 
dominant.  The  supercritical  mixing  and  combustion  of  LOX  and  methane  in  the  vicinity 
of  the  splitter  plate  has  been  carefully  analyzed.  The  study  confirms  that  the  flame  is 
stabilized  by  the  wake  recirculation  zone  immediately  downstream  of  the  splitter  plate. 
The  effect  of  real-fluid  thermodynamics  on  the  cryogen  flame  evolution  has  been 
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quantified.  The  strong  density  stratification  between  the  LOX  stream  and  hot  product 
exerts  a  substantial  influence  on  the  large-scale  flow  motions.  It  amplifies  the  axial 
turbulent  fluctuation  but  damps  the  radial  one,  and  thus,  leads  to  the  strong  anisotropy  of 
turbulence  near  the  density  interface. 
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